File indexing completed on 2024-04-06 12:02:18
0001
0002
0003
0004
0005
0006
0007 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0008 #include "CondFormats/JetMETObjects/interface/JetCorrectorParametersHelper.h"
0009 #include <iostream>
0010 #include <iomanip>
0011 #include <fstream>
0012 #include <sstream>
0013 #include <cstdlib>
0014 #include <algorithm>
0015 #include <cmath>
0016 #include <iterator>
0017
0018
0019
0020
0021
0022 JetCorrectorParameters::Definitions::Definitions(const std::vector<std::string>& fBinVar,
0023 const std::vector<std::string>& fParVar,
0024 const std::string& fFormula,
0025 bool fIsResponse) {
0026 for (unsigned i = 0; i < fBinVar.size(); i++)
0027 mBinVar.push_back(fBinVar[i]);
0028 for (unsigned i = 0; i < fParVar.size(); i++)
0029 mParVar.push_back(fParVar[i]);
0030 mFormula = fFormula;
0031 mIsResponse = fIsResponse;
0032 mLevel = "";
0033 }
0034
0035
0036
0037
0038 JetCorrectorParameters::Definitions::Definitions(const std::string& fLine) {
0039 std::vector<std::string> tokens = getTokens(fLine);
0040 if (!tokens.empty()) {
0041 if (tokens.size() < 6) {
0042 std::stringstream sserr;
0043 sserr << "(line " << fLine << "): less than 6 expected tokens:" << tokens.size();
0044 handleError("JetCorrectorParameters::Definitions", sserr.str());
0045 }
0046 unsigned nvar = getUnsigned(tokens[0]);
0047 unsigned npar = getUnsigned(tokens[nvar + 1]);
0048 for (unsigned i = 0; i < nvar; i++)
0049 mBinVar.push_back(tokens[i + 1]);
0050 for (unsigned i = 0; i < npar; i++)
0051 mParVar.push_back(tokens[nvar + 2 + i]);
0052 mFormula = tokens[npar + nvar + 2];
0053 std::string ss = tokens[npar + nvar + 3];
0054 if (ss == "Response")
0055 mIsResponse = true;
0056 else if (ss == "Correction")
0057 mIsResponse = false;
0058 else if (ss == "Resolution")
0059 mIsResponse = false;
0060 else if (ss.find("PAR") == 0)
0061 mIsResponse = false;
0062 else {
0063 std::stringstream sserr;
0064 sserr << "unknown option (" << ss << ")";
0065 handleError("JetCorrectorParameters::Definitions", sserr.str());
0066 }
0067 mLevel = tokens[npar + nvar + 4];
0068 }
0069 }
0070
0071
0072
0073
0074 JetCorrectorParameters::Record::Record(const std::string& fLine, unsigned fNvar) : mMin(0), mMax(0) {
0075 mNvar = fNvar;
0076
0077 std::vector<std::string> tokens = getTokens(fLine);
0078 if (!tokens.empty()) {
0079 if (tokens.size() < 3) {
0080 std::stringstream sserr;
0081 sserr << "(line " << fLine << "): "
0082 << "three tokens expected, " << tokens.size() << " provided.";
0083 handleError("JetCorrectorParameters::Record", sserr.str());
0084 }
0085 for (unsigned i = 0; i < mNvar; i++) {
0086 mMin.push_back(getFloat(tokens[i * 2]));
0087 mMax.push_back(getFloat(tokens[i * 2 + 1]));
0088 }
0089 unsigned nParam = getUnsigned(tokens[2 * mNvar]);
0090 if (nParam != tokens.size() - (2 * mNvar + 1)) {
0091 std::stringstream sserr;
0092 sserr << "(line " << fLine << "): " << tokens.size() - (2 * mNvar + 1) << " parameters, but nParam=" << nParam
0093 << ".";
0094 handleError("JetCorrectorParameters::Record", sserr.str());
0095 }
0096 for (unsigned i = (2 * mNvar + 1); i < tokens.size(); ++i)
0097 mParameters.push_back(getFloat(tokens[i]));
0098 }
0099 }
0100 std::ostream& operator<<(std::ostream& out, const JetCorrectorParameters::Record& fBin) {
0101 for (unsigned j = 0; j < fBin.nVar(); j++)
0102 out << fBin.xMin(j) << " " << fBin.xMax(j) << " ";
0103 out << fBin.nParameters() << " ";
0104 for (unsigned j = 0; j < fBin.nParameters(); j++)
0105 out << fBin.parameter(j) << " ";
0106 out << std::endl;
0107 return out;
0108 }
0109
0110
0111
0112
0113 JetCorrectorParameters::JetCorrectorParameters(const std::string& fFile, const std::string& fSection) {
0114 std::ifstream input(fFile.c_str());
0115 std::string currentSection = "";
0116 std::string line;
0117 std::string currentDefinitions = "";
0118 while (std::getline(input, line)) {
0119 if (line.empty())
0120 continue;
0121 std::string section = getSection(line);
0122 std::string tmp = getDefinitions(line);
0123 if (!section.empty() && tmp.empty()) {
0124 currentSection = section;
0125 continue;
0126 }
0127 if (currentSection == fSection) {
0128 if (!tmp.empty()) {
0129 currentDefinitions = tmp;
0130 continue;
0131 }
0132 Definitions definitions(currentDefinitions);
0133 if (!(definitions.nBinVar() == 0 && definitions.formula().empty()))
0134 mDefinitions = definitions;
0135 Record record(line, mDefinitions.nBinVar());
0136 bool check(true);
0137 for (unsigned i = 0; i < mDefinitions.nBinVar(); ++i)
0138 if (record.xMin(i) == 0 && record.xMax(i) == 0)
0139 check = false;
0140 if (record.nParameters() == 0)
0141 check = false;
0142 if (check)
0143 mRecords.push_back(record);
0144 }
0145 }
0146 if (currentDefinitions.empty())
0147 handleError("JetCorrectorParameters", "No definitions found!!!");
0148 if (mRecords.empty() && currentSection.empty())
0149 mRecords.push_back(Record());
0150 if (mRecords.empty() && !currentSection.empty()) {
0151 std::stringstream sserr;
0152 sserr << "the requested section " << fSection << " doesn't exist!";
0153 handleError("JetCorrectorParameters", sserr.str());
0154 }
0155 std::sort(mRecords.begin(), mRecords.end());
0156 valid_ = true;
0157
0158 if (mDefinitions.nBinVar() <= MAX_SIZE_DIMENSIONALITY) {
0159 init();
0160 } else {
0161 std::stringstream sserr;
0162 sserr << "since the binned dimensionality is greater than " << MAX_SIZE_DIMENSIONALITY
0163 << " the SimpleJetCorrector will default to using the legacy binIndex function!";
0164 handleError("JetCorrectorParameters", sserr.str());
0165 }
0166 }
0167
0168
0169
0170 void JetCorrectorParameters::init() {
0171 std::sort(mRecords.begin(), mRecords.end());
0172 helper = std::make_shared<JetCorrectorParametersHelper>();
0173 helper->init(mDefinitions, mRecords);
0174 }
0175
0176
0177
0178 int JetCorrectorParameters::binIndex(const std::vector<float>& fX) const {
0179 int result = -1;
0180 unsigned N = mDefinitions.nBinVar();
0181 if (N != fX.size()) {
0182 std::stringstream sserr;
0183 sserr << "# bin variables " << N << " doesn't correspont to requested #: " << fX.size();
0184 handleError("JetCorrectorParameters", sserr.str());
0185 }
0186 unsigned tmp;
0187 for (unsigned i = 0; i < size(); ++i) {
0188 tmp = 0;
0189 for (unsigned j = 0; j < N; j++)
0190 if (fX[j] >= record(i).xMin(j) && fX[j] < record(i).xMax(j))
0191 tmp += 1;
0192 if (tmp == N) {
0193 result = i;
0194 break;
0195 }
0196 }
0197 return result;
0198 }
0199
0200
0201
0202 int JetCorrectorParameters::binIndexN(const std::vector<float>& fX) const {
0203 if (helper->size() > 0) {
0204 return helper->binIndexN(fX, mRecords);
0205 } else {
0206 return -1;
0207 }
0208 }
0209
0210
0211
0212 int JetCorrectorParameters::neighbourBin(unsigned fIndex, unsigned fVar, bool fNext) const {
0213 int result = -1;
0214 unsigned N = mDefinitions.nBinVar();
0215 if (fVar >= N) {
0216 std::stringstream sserr;
0217 sserr << "# of bin variables " << N << " doesn't correspond to requested #: " << fVar;
0218 handleError("JetCorrectorParameters", sserr.str());
0219 }
0220 unsigned tmp;
0221 for (unsigned i = 0; i < size(); ++i) {
0222 tmp = 0;
0223 for (unsigned j = 0; j < fVar; j++)
0224 if (fabs(record(i).xMin(j) - record(fIndex).xMin(j)) < 0.0001)
0225 tmp += 1;
0226 for (unsigned j = fVar + 1; j < N; j++)
0227 if (fabs(record(i).xMin(j) - record(fIndex).xMin(j)) < 0.0001)
0228 tmp += 1;
0229 if (tmp < N - 1)
0230 continue;
0231 if (tmp == N - 1) {
0232 if (fNext)
0233 if (fabs(record(i).xMin(fVar) - record(fIndex).xMax(fVar)) < 0.0001)
0234 tmp += 1;
0235 if (!fNext)
0236 if (fabs(record(i).xMax(fVar) - record(fIndex).xMin(fVar)) < 0.0001)
0237 tmp += 1;
0238 }
0239 if (tmp == N) {
0240 result = i;
0241 break;
0242 }
0243 }
0244 return result;
0245 }
0246
0247
0248
0249 unsigned JetCorrectorParameters::size(unsigned fVar) const {
0250 if (fVar >= mDefinitions.nBinVar()) {
0251 std::stringstream sserr;
0252 sserr << "requested bin variable index " << fVar << " is greater than number of variables "
0253 << mDefinitions.nBinVar();
0254 handleError("JetCorrectorParameters", sserr.str());
0255 }
0256 unsigned result = 0;
0257 float tmpMin(-9999), tmpMax(-9999);
0258 for (unsigned i = 0; i < size(); ++i)
0259 if (record(i).xMin(fVar) > tmpMin && record(i).xMax(fVar) > tmpMax) {
0260 result++;
0261 tmpMin = record(i).xMin(fVar);
0262 tmpMax = record(i).xMax(fVar);
0263 }
0264 return result;
0265 }
0266
0267
0268
0269 std::vector<float> JetCorrectorParameters::binCenters(unsigned fVar) const {
0270 std::vector<float> result;
0271 for (unsigned i = 0; i < size(); ++i)
0272 result.push_back(record(i).xMiddle(fVar));
0273 return result;
0274 }
0275
0276
0277
0278 void JetCorrectorParameters::printScreen() const {
0279 std::cout << "--------------------------------------------" << std::endl;
0280 std::cout << "//////// PARAMETERS: //////////////////////" << std::endl;
0281 std::cout << "--------------------------------------------" << std::endl;
0282 std::cout << "Number of binning variables: " << definitions().nBinVar() << std::endl;
0283 std::cout << "Names of binning variables: ";
0284 for (unsigned i = 0; i < definitions().nBinVar(); i++)
0285 std::cout << definitions().binVar(i) << " ";
0286 std::cout << std::endl;
0287 std::cout << "--------------------------------------------" << std::endl;
0288 std::cout << "Number of parameter variables: " << definitions().nParVar() << std::endl;
0289 std::cout << "Names of parameter variables: ";
0290 for (unsigned i = 0; i < definitions().nParVar(); i++)
0291 std::cout << definitions().parVar(i) << " ";
0292 std::cout << std::endl;
0293 std::cout << "--------------------------------------------" << std::endl;
0294 std::cout << "Parametrization Formula: " << definitions().formula() << std::endl;
0295 if (definitions().isResponse())
0296 std::cout << "Type (Response or Correction): "
0297 << "Response" << std::endl;
0298 else
0299 std::cout << "Type (Response or Correction): "
0300 << "Correction" << std::endl;
0301 std::cout << "Correction Level: " << definitions().level() << std::endl;
0302 std::cout << "--------------------------------------------" << std::endl;
0303 std::cout << "------- Bin contents -----------------------" << std::endl;
0304 for (unsigned i = 0; i < size(); i++) {
0305 for (unsigned j = 0; j < definitions().nBinVar(); j++)
0306 std::cout << record(i).xMin(j) << " " << record(i).xMax(j) << " ";
0307 std::cout << record(i).nParameters() << " ";
0308 for (unsigned j = 0; j < record(i).nParameters(); j++)
0309 std::cout << record(i).parameter(j) << " ";
0310 std::cout << std::endl;
0311 }
0312 }
0313
0314
0315
0316 void JetCorrectorParameters::printFile(const std::string& fFileName) const {
0317 std::ofstream txtFile;
0318 txtFile.open(fFileName.c_str());
0319 txtFile.setf(std::ios::right);
0320 txtFile << "{" << definitions().nBinVar() << std::setw(15);
0321 for (unsigned i = 0; i < definitions().nBinVar(); i++)
0322 txtFile << definitions().binVar(i) << std::setw(15);
0323 txtFile << definitions().nParVar() << std::setw(15);
0324 for (unsigned i = 0; i < definitions().nParVar(); i++)
0325 txtFile << definitions().parVar(i) << std::setw(15);
0326 txtFile << std::setw(definitions().formula().size() + 15) << definitions().formula() << std::setw(15);
0327 if (definitions().isResponse())
0328 txtFile << "Response" << std::setw(15);
0329 else
0330 txtFile << "Correction" << std::setw(15);
0331 txtFile << definitions().level() << "}"
0332 << "\n";
0333 for (unsigned i = 0; i < size(); i++) {
0334 for (unsigned j = 0; j < definitions().nBinVar(); j++)
0335 txtFile << record(i).xMin(j) << std::setw(15) << record(i).xMax(j) << std::setw(15);
0336 txtFile << record(i).nParameters() << std::setw(15);
0337 for (unsigned j = 0; j < record(i).nParameters(); j++)
0338 txtFile << record(i).parameter(j) << std::setw(15);
0339 txtFile << "\n";
0340 }
0341 txtFile.close();
0342 }
0343
0344 namespace {
0345 const std::vector<std::string> labels_ = {
0346 "L1Offset",
0347 "L2Relative",
0348 "L3Absolute",
0349 "L4EMF",
0350 "L5Flavor",
0351 "L6UE",
0352 "L7Parton",
0353 "L1JPTOffset",
0354 "L2L3Residual",
0355 "Uncertainty",
0356 "L1FastJet",
0357 "UncertaintyAbsolute",
0358 "UncertaintyHighPtExtra",
0359 "UncertaintySinglePionECAL",
0360 "UncertaintyFlavor",
0361 "UncertaintyTime",
0362 "UncertaintyRelativeJEREC1",
0363 "UncertaintyRelativeJEREC2",
0364 "UncertaintyRelativeJERHF",
0365 "UncertaintyRelativeStatEC2",
0366 "UncertaintyRelativeStatHF",
0367 "UncertaintyRelativeFSR",
0368 "UncertaintyPileUpDataMC",
0369 "UncertaintyPileUpOOT",
0370 "UncertaintyPileUpPtBB",
0371 "UncertaintyPileUpBias",
0372 "UncertaintyPileUpJetRate",
0373 "UncertaintySinglePionHCAL",
0374 "UncertaintyRelativePtEC1",
0375 "UncertaintyRelativePtEC2",
0376 "UncertaintyRelativePtHF",
0377 "UncertaintyRelativeSample",
0378 "UncertaintyPileUpPtEC",
0379 "UncertaintyPileUpPtHF",
0380 "L1RC",
0381 "L1Residual",
0382 "UncertaintyAux3",
0383 "UncertaintyAux4",
0384 };
0385
0386 const std::vector<std::string> l5Flavors_ = {"L5Flavor_bJ",
0387 "L5Flavor_cJ",
0388 "L5Flavor_qJ",
0389 "L5Flavor_gJ",
0390 "L5Flavor_bT",
0391 "L5Flavor_cT",
0392 "L5Flavor_qT",
0393 "L5Flavor_gT"};
0394
0395 const std::vector<std::string> l7Partons_ = {"L7Parton_gJ",
0396 "L7Parton_qJ",
0397 "L7Parton_cJ",
0398 "L7Parton_bJ",
0399 "L7Parton_jJ",
0400 "L7Parton_qT",
0401 "L7Parton_cT",
0402 "L7Parton_bT",
0403 "L7Parton_tT"};
0404 }
0405 std::string JetCorrectorParametersCollection::findLabel(key_type k) {
0406 if (isL5(k))
0407 return findL5Flavor(k);
0408 else if (isL7(k))
0409 return findL7Parton(k);
0410 else
0411 return labels_[k];
0412 }
0413
0414 std::string JetCorrectorParametersCollection::findL5Flavor(key_type k) {
0415 if (k == L5Flavor)
0416 return labels_[L5Flavor];
0417 else
0418 return l5Flavors_[k / 100 - 1];
0419 }
0420
0421 std::string JetCorrectorParametersCollection::findL7Parton(key_type k) {
0422 if (k == L7Parton)
0423 return labels_[L7Parton];
0424 else
0425 return l7Partons_[k / 1000 - 1];
0426 }
0427
0428 void JetCorrectorParametersCollection::getSections(std::string inputFile, std::vector<std::string>& outputs) {
0429 outputs.clear();
0430 std::ifstream input(inputFile.c_str());
0431 while (!input.eof()) {
0432 char buff[10000];
0433 input.getline(buff, 10000);
0434 std::string in(buff);
0435 if (in[0] == '[') {
0436 std::string tok = getSection(in);
0437 if (!tok.empty()) {
0438 outputs.push_back(tok);
0439 }
0440 }
0441 }
0442 std::cout << "Found these sections for file: " << std::endl;
0443 copy(outputs.begin(), outputs.end(), std::ostream_iterator<std::string>(std::cout, "\n"));
0444 }
0445
0446
0447 void JetCorrectorParametersCollection::push_back(key_type i, value_type const& j, label_type const& flav) {
0448 std::cout << "i = " << i << std::endl;
0449 std::cout << "flav = " << flav << std::endl;
0450 if (isL5(i)) {
0451 std::cout << "This is L5, getL5Bin = " << getL5Bin(flav) << std::endl;
0452 correctionsL5_.push_back(pair_type(getL5Bin(flav), j));
0453 } else if (isL7(i)) {
0454 std::cout << "This is L7, getL7Bin = " << getL7Bin(flav) << std::endl;
0455 correctionsL7_.push_back(pair_type(getL7Bin(flav), j));
0456 } else if (flav.empty()) {
0457 corrections_.push_back(pair_type(i, j));
0458 } else {
0459 std::cout << "***** NOT ADDING " << flav << ", corresponding position in JetCorrectorParameters is not found."
0460 << std::endl;
0461 }
0462 }
0463
0464
0465
0466 JetCorrectorParameters const& JetCorrectorParametersCollection::operator[](key_type k) const {
0467 collection_type::const_iterator ibegin, iend, i;
0468 if (isL5(k)) {
0469 ibegin = correctionsL5_.begin();
0470 iend = correctionsL5_.end();
0471 i = ibegin;
0472 } else if (isL7(k)) {
0473 ibegin = correctionsL7_.begin();
0474 iend = correctionsL7_.end();
0475 i = ibegin;
0476 } else {
0477 ibegin = corrections_.begin();
0478 iend = corrections_.end();
0479 i = ibegin;
0480 }
0481 for (; i != iend; ++i) {
0482 if (k == i->first)
0483 return i->second;
0484 }
0485 throw cms::Exception("InvalidInput") << " cannot find key " << static_cast<int>(k)
0486 << " in the JEC payload, this usually means you have to change the global tag"
0487 << std::endl;
0488 }
0489
0490
0491
0492 void JetCorrectorParametersCollection::validKeys(std::vector<key_type>& keys) const {
0493 keys.clear();
0494 for (collection_type::const_iterator ibegin = corrections_.begin(), iend = corrections_.end(), i = ibegin; i != iend;
0495 ++i) {
0496 keys.push_back(i->first);
0497 }
0498 for (collection_type::const_iterator ibegin = correctionsL5_.begin(), iend = correctionsL5_.end(), i = ibegin;
0499 i != iend;
0500 ++i) {
0501 keys.push_back(i->first);
0502 }
0503 for (collection_type::const_iterator ibegin = correctionsL7_.begin(), iend = correctionsL7_.end(), i = ibegin;
0504 i != iend;
0505 ++i) {
0506 keys.push_back(i->first);
0507 }
0508 }
0509
0510
0511 JetCorrectorParametersCollection::key_type JetCorrectorParametersCollection::getL5Bin(std::string const& flav) {
0512 std::vector<std::string>::const_iterator found = find(l5Flavors_.begin(), l5Flavors_.end(), flav);
0513 if (found != l5Flavors_.end()) {
0514 return (found - l5Flavors_.begin() + 1) * 100;
0515 } else
0516 return L5Flavor;
0517 }
0518
0519 JetCorrectorParametersCollection::key_type JetCorrectorParametersCollection::getL7Bin(std::string const& flav) {
0520 std::vector<std::string>::const_iterator found = find(l7Partons_.begin(), l7Partons_.end(), flav);
0521 if (found != l7Partons_.end()) {
0522 return (found - l7Partons_.begin() + 1) * 1000;
0523 } else
0524 return L7Parton;
0525 }
0526
0527
0528 bool JetCorrectorParametersCollection::isL5(key_type k) { return k == L5Flavor || (k / 100 > 0 && k / 1000 == 0); }
0529
0530 bool JetCorrectorParametersCollection::isL7(key_type k) { return k == L7Parton || (k / 1000 > 0); }
0531
0532
0533 JetCorrectorParametersCollection::key_type JetCorrectorParametersCollection::findKey(std::string const& label) const {
0534
0535 std::vector<std::string>::const_iterator found1 = find(l5Flavors_.begin(), l5Flavors_.end(), label);
0536 if (found1 != l5Flavors_.end()) {
0537 return getL5Bin(label);
0538 }
0539
0540
0541 std::vector<std::string>::const_iterator found2 = find(l7Partons_.begin(), l7Partons_.end(), label);
0542 if (found2 != l7Partons_.end()) {
0543 return getL7Bin(label);
0544 }
0545
0546
0547 std::vector<std::string>::const_iterator found3 = find(labels_.begin(), labels_.end(), label);
0548 if (found3 != labels_.end()) {
0549 return static_cast<key_type>(found3 - labels_.begin());
0550 }
0551
0552
0553 throw cms::Exception("InvalidInput") << " Cannot find label " << label << std::endl;
0554 }
0555
0556
0557
0558
0559 #include "FWCore/Utilities/interface/typelookup.h"
0560
0561 TYPELOOKUP_DATA_REG(JetCorrectorParameters);
0562 TYPELOOKUP_DATA_REG(JetCorrectorParametersCollection);