Project CMSSW displayed by LXR CMSSW_15_0_X_2024-12-02-2300/​CondFormats/​JetMETObjects/​src/​JetResolutionObject.cc	Source navigation Diff markup Identifier search General search
	Version: CMSSW_15_0_X_2024-12-02-2300 CMSSW_15_0_X_2024-12-01-2300 CMSSW_15_0_X_2024-11-28-2300 CMSSW_15_0_X_2024-11-27-2300 CMSSW_15_0_X_2024-11-26-2300 Revert   Change

File indexing completed on 2024-04-06 12:02:18

 #ifndef STANDALONE
 #include <CondFormats/JetMETObjects/interface/JetResolutionObject.h>
 #include <CondFormats/JetMETObjects/interface/Utilities.h>
 #include <FWCore/Utilities/interface/EDMException.h>
 
 #else
 #include "JetResolutionObject.h"
 #include "Utilities.h"
 #include <exception>
 
 namespace edm {
   namespace errors {
     enum ErrorCode { NotFound = 8026, ConfigFileReadError = 7002, UnimplementedFeature = 8011, FileReadError = 8021 };
   };
 };  // namespace edm
 #endif
 
 #include <cmath>
 #include <iostream>
 #include <fstream>
 #include <iomanip>
 #include <algorithm>
 
 namespace JME {
 
   std::string getDefinitionLine(const std::string& line) {
     size_t first = line.find('{');
     size_t last = line.find('}');
 
     if (first != std::string::npos && last != std::string::npos && first < last)
       return std::string(line, first + 1, last - first - 1);
 
     return "";
   }
 
   void throwException(uint32_t code, const std::string& message) {
 #ifndef STANDALONE
     throw edm::Exception(static_cast<edm::errors::ErrorCodes>(code), message);
 #else
     std::stringstream error;
     error << message << " Error code: " << code;
     throw std::runtime_error(error.str());
 
 #endif
   }
 
   const bimap<Binning, std::string> JetParameters::binning_to_string = {{Binning::JetPt, "JetPt"},
                                                                         {Binning::JetEta, "JetEta"},
                                                                         {Binning::JetAbsEta, "JetAbsEta"},
                                                                         {Binning::JetE, "JetE"},
                                                                         {Binning::JetArea, "JetArea"},
                                                                         {Binning::Mu, "Mu"},
                                                                         {Binning::Rho, "Rho"},
                                                                         {Binning::NPV, "NPV"}};
 
   JetParameters::JetParameters(JetParameters&& rhs) { m_values = std::move(rhs.m_values); }
 
   JetParameters::JetParameters(std::initializer_list<typename value_type::value_type> init) {
     for (auto& i : init) {
       set(i.first, i.second);
     }
   }
 
   JetParameters& JetParameters::setJetPt(float pt) {
     m_values[Binning::JetPt] = pt;
     return *this;
   }
 
   JetParameters& JetParameters::setJetEta(float eta) {
     m_values[Binning::JetEta] = eta;
     m_values[Binning::JetAbsEta] = fabs(eta);
     return *this;
   }
 
   JetParameters& JetParameters::setJetE(float e) {
     m_values[Binning::JetE] = e;
     return *this;
   }
 
   JetParameters& JetParameters::setJetArea(float area) {
     m_values[Binning::JetArea] = area;
     return *this;
   }
 
   JetParameters& JetParameters::setMu(float mu) {
     m_values[Binning::Mu] = mu;
     return *this;
   }
 
   JetParameters& JetParameters::setNPV(float npv) {
     m_values[Binning::NPV] = npv;
     return *this;
   }
 
   JetParameters& JetParameters::setRho(float rho) {
     m_values[Binning::Rho] = rho;
     return *this;
   }
 
   JetParameters& JetParameters::set(const Binning& bin, float value) {
     m_values.emplace(bin, value);
 
     // Special case for eta
     if (bin == Binning::JetEta) {
       m_values.emplace(Binning::JetAbsEta, fabs(value));
     }
 
     return *this;
   }
 
   JetParameters& JetParameters::set(const typename value_type::value_type& value) {
     set(value.first, value.second);
     return *this;
   }
 
   std::vector<float> JetParameters::createVector(const std::vector<Binning>& binning) const {
     std::vector<float> values;
 
     for (const auto& bin : binning) {
       const auto& it = m_values.find(bin);
       if (it == m_values.cend()) {
         throwException(edm::errors::NotFound,
                        "JER parametrisation depends on '" + JetParameters::binning_to_string.left.at(bin) +
                            "' but no value for this parameter has been specified. Please call the appropriate 'set' "
                            "function of the JME::JetParameters object");
       }
 
       values.push_back(it->second);
     }
 
     return values;
   }
 
   std::vector<float> JetParameters::createVector(const std::vector<std::string>& binname) const {
     std::vector<float> values;
 
     for (const auto& name : binname) {
       Binning bi = binning_to_string.right.find(name)->second;
 
       const auto& it = m_values.find(bi);
       if (it == m_values.cend()) {
         edm::LogPrint("JPM") << "Bin name " << name << " not found!";
         throwException(edm::errors::NotFound,
                        "JER parametrisation depends on '" + JetParameters::binning_to_string.left.at(bi) +
                            "' but no value for this parameter has been specified. Please call the appropriate 'set' "
                            "function of the JME::JetParameters object");
       }
 
       values.push_back(it->second);
     }
 
     return values;
   }
 
   JetResolutionObject::Definition::Definition(const std::string& definition) {
     std::vector<std::string> tokens = getTokens(definition);
 
     // We need at least 3 tokens
     if (tokens.size() < 3) {
       throwException(edm::errors::ConfigFileReadError,
                      "Definition line needs at least three tokens. Please check file format.");
     }
 
     size_t n_bins = std::stoul(tokens[0]);
 
     if (tokens.size() < (n_bins + 2)) {
       throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check.");
     }
 
     for (size_t i = 0; i < n_bins; i++) {
       m_bins_name.push_back(tokens[i + 1]);
     }
 
     size_t n_variables = std::stoul(tokens[n_bins + 1]);
 
     if (tokens.size() < (1 + n_bins + 1 + n_variables + 1)) {
       throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check.");
     }
 
     for (size_t i = 0; i < n_variables; i++) {
       m_variables_name.push_back(tokens[n_bins + 2 + i]);
     }
 
     m_formula_str = tokens[n_bins + n_variables + 2];
 
     std::string formula_str_lower = m_formula_str;
     std::transform(formula_str_lower.begin(), formula_str_lower.end(), formula_str_lower.begin(), ::tolower);
 
     if (formula_str_lower == "none") {
       m_formula_str = "";
 
       if ((tokens.size() > n_bins + n_variables + 3) && (std::atoi(tokens[n_bins + n_variables + 3].c_str()))) {
         size_t n_parameters = std::stoul(tokens[n_bins + n_variables + 3]);
 
         if (tokens.size() < (1 + n_bins + 1 + n_variables + 1 + 1 + n_parameters)) {
           throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check.");
         }
 
         for (size_t i = 0; i < n_parameters; i++) {
           m_formula_str += tokens[n_bins + n_variables + 4 + i] + " ";
         }
       }
     }
 
     init();
   }
 
   void JetResolutionObject::Definition::init() {
     if (!m_formula_str.empty()) {
       if (m_formula_str.find(' ') == std::string::npos)
 #ifndef STANDALONE
         m_formula = std::make_shared<reco::FormulaEvaluator>(m_formula_str);
 #else
         m_formula = std::make_shared<TFormula>("jet_resolution_formula", m_formula_str.c_str());
 #endif
       else
         m_parameters_name = getTokens(m_formula_str);
     }
     for (const auto& bin : m_bins_name) {
       const auto& b = JetParameters::binning_to_string.right.find(bin);
       if (b == JetParameters::binning_to_string.right.cend()) {
         throwException(edm::errors::UnimplementedFeature, "Bin name not supported: '" + bin + "'");
       }
       m_bins.push_back(b->second);
     }
 
     for (const auto& v : m_variables_name) {
       const auto& var = JetParameters::binning_to_string.right.find(v);
       if (var == JetParameters::binning_to_string.right.cend()) {
         throwException(edm::errors::UnimplementedFeature, "Variable name not supported: '" + v + "'");
       }
       m_variables.push_back(var->second);
     }
   }
 
   JetResolutionObject::Record::Record(const std::string& line, const Definition& def) {
     std::vector<std::string> tokens = getTokens(line);
 
     if (tokens.size() < (def.nBins() * 2 + def.nVariables() * 2 + 1)) {
       throwException(edm::errors::ConfigFileReadError, "Invalid record. Please check file format. Record: " + line);
     }
 
     size_t pos = 0;
 
     for (size_t i = 0; i < def.nBins(); i++) {
       Range r(std::stof(tokens[pos]), std::stof(tokens[pos + 1]));
       pos += 2;
       m_bins_range.push_back(r);
     }
 
     size_t n_parameters = std::stoul(tokens[pos++]);
 
     if (tokens.size() < (def.nBins() * 2 + def.nVariables() * 2 + 1 + (n_parameters - def.nVariables() * 2))) {
       throwException(edm::errors::ConfigFileReadError, "Invalid record. Please check file format. Record: " + line);
     }
 
     for (size_t i = 0; i < def.nVariables(); i++) {
       Range r(std::stof(tokens[pos]), std::stof(tokens[pos + 1]));
       pos += 2;
       m_variables_range.push_back(r);
       n_parameters -= 2;
     }
 
     for (size_t i = 0; i < n_parameters; i++) {
       m_parameters_values.push_back(std::stof(tokens[pos++]));
     }
   }
 
   JetResolutionObject::JetResolutionObject(const std::string& filename) {
     // Parse file
     std::ifstream f(filename);
 
     if (!f.good()) {
       throwException(edm::errors::FileReadError, "Can't read input file '" + filename + "'");
     }
 
     for (std::string line; std::getline(f, line);) {
       if ((line.empty()) || (line[0] == '#'))
         continue;
 
       std::string definition = getDefinitionLine(line);
 
       if (!definition.empty()) {
         m_definition = Definition(definition);
       } else {
         m_records.push_back(Record(line, m_definition));
       }
     }
 
     m_valid = true;
   }
 
   JetResolutionObject::JetResolutionObject(const JetResolutionObject& object) {
     m_definition = object.m_definition;
     m_records = object.m_records;
     m_valid = object.m_valid;
 
     m_definition.init();
   }
 
   JetResolutionObject::JetResolutionObject() {
     // Empty
   }
 
   void JetResolutionObject::dump() const {
     std::cout << "Definition: " << std::endl;
     std::cout << "    Number of binning variables: " << m_definition.nBins() << std::endl;
     std::cout << "        ";
     for (const auto& bin : m_definition.getBinsName()) {
       std::cout << bin << ", ";
     }
     std::cout << std::endl;
     std::cout << "    Number of variables: " << m_definition.nVariables() << std::endl;
     std::cout << "        ";
     for (const auto& bin : m_definition.getVariablesName()) {
       std::cout << bin << ", ";
     }
     std::cout << std::endl;
     std::cout << "    Formula: " << m_definition.getFormulaString() << std::endl;
 
     std::cout << std::endl << "Bin contents" << std::endl;
 
     for (const auto& record : m_records) {
       std::cout << "    Bins" << std::endl;
       size_t index = 0;
       for (const auto& bin : record.getBinsRange()) {
         std::cout << "        " << m_definition.getBinName(index) << " [" << bin.min << " - " << bin.max << "]"
                   << std::endl;
         index++;
       }
 
       std::cout << "    Variables" << std::endl;
       index = 0;
       for (const auto& r : record.getVariablesRange()) {
         std::cout << "        " << m_definition.getVariableName(index) << " [" << r.min << " - " << r.max << "] "
                   << std::endl;
         index++;
       }
 
       std::cout << "    Parameters" << std::endl;
       index = 0;
       for (const auto& par : record.getParametersValues()) {
         std::cout << "        Parameter #" << index << " = " << par << std::endl;
         index++;
       }
     }
   }
 
   void JetResolutionObject::saveToFile(const std::string& file) const {
     std::ofstream fout(file);
     fout.setf(std::ios::right);
 
     // Definition
     fout << "{" << m_definition.nBins();
 
     for (auto& bin : m_definition.getBinsName())
       fout << "    " << bin;
 
     fout << "    " << m_definition.nVariables();
 
     for (auto& var : m_definition.getVariablesName())
       fout << "    " << var;
 
     fout << "    " << (m_definition.getFormulaString().empty() ? "None" : m_definition.getFormulaString())
          << "    Resolution}" << std::endl;
 
     // Records
     for (auto& record : m_records) {
       for (auto& r : record.getBinsRange()) {
         fout << std::left << std::setw(15) << r.min << std::setw(15) << r.max << std::setw(15);
       }
       fout << (record.nVariables() * 2 + record.nParameters()) << std::setw(15);
 
       for (auto& r : record.getVariablesRange()) {
         fout << r.min << std::setw(15) << r.max << std::setw(15);
       }
 
       for (auto& p : record.getParametersValues()) {
         fout << p << std::setw(15);
       }
 
       fout << std::endl << std::setw(0);
     }
   }
 
   const JetResolutionObject::Record* JetResolutionObject::getRecord(const JetParameters& bins_parameters) const {
     // Find record for bins
     if (!m_valid)
       return nullptr;
 
     // Create vector of bins value. Throw if some values are missing
     std::vector<float> bins = bins_parameters.createVector(m_definition.getBinsName());
 
     // Iterate over all records, and find the one for which all bins are valid
     const Record* good_record = nullptr;
     for (const auto& record : m_records) {
       // Iterate over bins
       size_t valid_bins = 0;
       size_t current_bin = 0;
       for (const auto& bin : record.getBinsRange()) {
         if (bin.is_inside(bins[current_bin])) {
           valid_bins++;
         }
 
         current_bin++;
       }
 
       if (valid_bins == m_definition.nBins()) {
         good_record = &record;
         break;
       }
     }
 
     return good_record;
   }
 
   float JetResolutionObject::evaluateFormula(const JetResolutionObject::Record& record,
                                              const JetParameters& variables_parameters) const {
     if (!m_valid)
       return 1;
 
 #ifndef STANDALONE
     reco::FormulaEvaluator formula = reco::FormulaEvaluator(m_definition.getFormulaString());
 #else
     // Set parameters
     auto const* pFormula = m_definition.getFormula();
     if (!pFormula)
       return 1;
     auto formula = *pFormula;
 #endif
     // Create vector of variables value. Throw if some values are missing
     std::vector<float> variables = variables_parameters.createVector(m_definition.getVariablesName());
 
     double variables_[4] = {0};
     for (size_t index = 0; index < m_definition.nVariables(); index++) {
       variables_[index] =
           clip(variables[index], record.getVariablesRange()[index].min, record.getVariablesRange()[index].max);
     }
 
     const std::vector<float>& parameters = record.getParametersValues();
 
 #ifndef STANDALONE
     //ArrayAdaptor only takes doubles
     std::vector<double> parametersD(parameters.begin(), parameters.end());
     return formula.evaluate(reco::formula::ArrayAdaptor(variables_, m_definition.nVariables()),
                             reco::formula::ArrayAdaptor(parametersD.data(), parametersD.size()));
 #else
     for (size_t index = 0; index < parameters.size(); index++) {
       formula.SetParameter(index, parameters[index]);
     }
 
     return formula.EvalPar(variables_);
 #endif
   }
 }  // namespace JME
 
 #ifndef STANDALONE
 #include "FWCore/Utilities/interface/typelookup.h"
 TYPELOOKUP_DATA_REG(JME::JetResolutionObject);
 #endif

This page was automatically generated by the 2.2.1 LXR engine. The LXR team