File indexing completed on 2023-03-17 10:47:07
0001
0002
0003
0004
0005
0006
0007 #include "CondFormats/JetMETObjects/interface/JetCorrectorParametersHelper.h"
0008 #include <sstream>
0009 #include <cstdlib>
0010 #include <algorithm>
0011 #include <cmath>
0012 #include <iterator>
0013
0014
0015
0016
0017
0018 void JetCorrectorParametersHelper::initTransientMaps() {
0019 mIndexMap.clear();
0020 mMap.clear();
0021 mBinBoundaries.clear();
0022 mBinBoundaries.assign(JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY, std::vector<float>(0, 0));
0023 }
0024 void JetCorrectorParametersHelper::init(const JetCorrectorParameters::Definitions& mDefinitionsLocal,
0025 const std::vector<JetCorrectorParameters::Record>& mRecordsLocal) {
0026 SIZE = mDefinitionsLocal.nBinVar();
0027 if (SIZE > JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY) {
0028 std::stringstream sserr;
0029 sserr << "The number of binned variables requested (" << SIZE << ") is greater than the number allowed ("
0030 << JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY << ")";
0031 handleError("JetCorrectorParametersHelper", sserr.str());
0032 }
0033
0034 initTransientMaps();
0035 size_t start = 0, end = 0;
0036 size_t nRec = mRecordsLocal.size();
0037 size_t indexMapSize = 0, tmpIndexMapSize = 0;
0038 for (unsigned i = 0; i < nRec; ++i) {
0039 for (unsigned j = 0; j < SIZE; j++) {
0040 if (j < SIZE - 1 && std::find(mBinBoundaries[j].begin(), mBinBoundaries[j].end(), mRecordsLocal[i].xMin(j)) ==
0041 mBinBoundaries[j].end())
0042 mBinBoundaries[j].push_back(mRecordsLocal[i].xMin(j));
0043 else if (j == SIZE - 1) {
0044 if (i == 0)
0045 mBinBoundaries[j].reserve(mRecordsLocal.size());
0046
0047 mBinBoundaries[j].push_back(mRecordsLocal[i].xMin(j));
0048
0049 if (SIZE > 1 && (i == nRec - 1 || mRecordsLocal[i].xMin(j - 1) != mRecordsLocal[i + 1].xMin(j - 1))) {
0050 end = i;
0051
0052 mMap.emplace(gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY - 1>(
0053 [&](size_t k) { return (k < SIZE - 1) ? mRecordsLocal[i].xMin(k) : -9999; }),
0054 std::make_pair(start, end));
0055 start = i + 1;
0056 }
0057 }
0058 }
0059 indexMapSize = mIndexMap.size();
0060 tuple_type tmpTuple = gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY>(
0061 [&](size_t k) { return (k < SIZE) ? mRecordsLocal[i].xMin(k) : -9999; });
0062 mIndexMap.emplace(tmpTuple, i);
0063 tmpIndexMapSize = mIndexMap.size();
0064 if (indexMapSize == tmpIndexMapSize) {
0065 size_t existing_index = mIndexMap.find(tmpTuple)->second;
0066 std::stringstream sserr;
0067 sserr << "Duplicate binning in record found (existing index,current index)=(" << existing_index << "," << i << ")"
0068 << std::endl
0069 << "\tBins(lower bounds)=" << tmpTuple;
0070 handleError("JetCorrectorParametersHelper", sserr.str());
0071 }
0072 }
0073 if (mBinBoundaries[SIZE - 1].size() != nRec) {
0074 std::stringstream sserr;
0075 sserr << "Did not find all bin boundaries for dimension " << SIZE - 1 << "!!!" << std::endl
0076 << "Found " << mBinBoundaries[SIZE - 1].size() << " out of " << nRec << " records";
0077 handleError("JetCorrectorParametersHelper", sserr.str());
0078 }
0079 indexMapSize = mIndexMap.size();
0080 if (indexMapSize != nRec) {
0081 handleError("JetCorrectorParametersHelper",
0082 "The mapping of bin lower bounds to indices does not contain all possible entries!!!");
0083 }
0084
0085
0086 if (SIZE > 2)
0087 checkMiddleBinUniformity(mRecordsLocal);
0088 }
0089
0090 void JetCorrectorParametersHelper::checkMiddleBinUniformity(
0091 const std::vector<JetCorrectorParameters::Record>& mRecords) const {
0092 unsigned N = SIZE - 2;
0093 size_t nRec = mRecords.size();
0094 std::vector<int> fN(N, -1);
0095
0096
0097 for (unsigned idim = 1; idim <= N; idim++) {
0098 int nBoundaryPassed = 0;
0099 if (fN[idim - 1] == -1)
0100 fN[idim - 1] = mBinBoundaries[idim].size();
0101
0102 for (unsigned iRecord = 0; iRecord < nRec; iRecord++) {
0103 if (mRecords[iRecord].xMin(idim) != mBinBoundaries[idim][nBoundaryPassed % fN[idim - 1]]) {
0104 throw cms::Exception("MissingRecord")
0105 << "found a missing record in binned dimension " << idim << " after record " << iRecord << std::endl
0106 << "\tthe bin lower bound should have been " << mBinBoundaries[idim][nBoundaryPassed % fN[idim - 1]]
0107 << ", but was instead " << mRecords[iRecord].xMin(idim) << std::endl
0108 << "\tall binned dimensions, besides the last one, must have uniform binning." << std::endl
0109 << mRecords[iRecord - 1] << std::endl
0110 << mRecords[iRecord] << std::endl;
0111 } else if ((iRecord == nRec - 1 || mRecords[iRecord].xMin(idim) != mRecords[iRecord + 1].xMin(idim))) {
0112 nBoundaryPassed++;
0113 } else {
0114 continue;
0115 }
0116 }
0117 }
0118 }
0119
0120
0121
0122
0123
0124 void JetCorrectorParametersHelper::binIndexChecks(unsigned N, const std::vector<float>& fX) const {
0125 if (N != fX.size()) {
0126 std::stringstream sserr;
0127 sserr << "The number of binned variables, " << N << ", doesn't correspond to the number requested, " << fX.size();
0128 handleError("JetCorrectorParametersHelper", sserr.str());
0129 }
0130 }
0131 bool JetCorrectorParametersHelper::binBoundChecks(unsigned dim,
0132 const float& value,
0133 const float& min,
0134 const float& max) const {
0135 if (value < min || value >= max)
0136 return false;
0137 else
0138 return true;
0139 }
0140
0141
0142
0143
0144 int JetCorrectorParametersHelper::binIndexN(const std::vector<float>& fX,
0145 const std::vector<JetCorrectorParameters::Record>& mRecords) const {
0146 unsigned Nm1 = SIZE - 1;
0147 binIndexChecks(SIZE, fX);
0148
0149
0150 std::vector<float> fN(SIZE, -1);
0151 std::vector<float>::const_iterator tmpIt;
0152
0153
0154 for (unsigned idim = 0; idim == 0 || idim < fX.size() - 1; idim++) {
0155 if (!binBoundChecks(idim, fX[idim], *mBinBoundaries[idim].begin(), mRecords[size() - 1].xMax(idim)))
0156 return -1;
0157 tmpIt = std::lower_bound(mBinBoundaries[idim].begin(), mBinBoundaries[idim].end(), fX[idim]);
0158
0159
0160 if (tmpIt == mBinBoundaries[idim].end())
0161 tmpIt = mBinBoundaries[idim].begin() + mBinBoundaries[idim].size() - 1;
0162 else if (*tmpIt != fX[idim])
0163 tmpIt -= 1;
0164 fN[idim] = *tmpIt;
0165 }
0166
0167
0168 std::pair<size_t, size_t> indexBounds;
0169 if (SIZE > 1) {
0170 tuple_type_Nm1 to_find_Nm1 = gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY - 1>(
0171 [&](size_t i) { return (i < Nm1) ? fN[i] : -9999; });
0172 if (mMap.find(to_find_Nm1) != mMap.end())
0173 indexBounds = mMap.at(to_find_Nm1);
0174 else {
0175 std::stringstream sserr;
0176 sserr << "couldn't find the index boundaries for dimension " << Nm1 << std::endl
0177 << "looking for last bin with N-1 values of " << to_find_Nm1 << std::endl;
0178 handleError("JetCorrectorParametersHelper", sserr.str());
0179 return -1;
0180 }
0181
0182
0183 if (!binBoundChecks(Nm1, fX[Nm1], mRecords[indexBounds.first].xMin(Nm1), mRecords[indexBounds.second].xMax(Nm1)))
0184 return -1;
0185 tmpIt = std::lower_bound(
0186 mBinBoundaries[Nm1].begin() + indexBounds.first, mBinBoundaries[Nm1].begin() + indexBounds.second, fX[Nm1]);
0187 if (*tmpIt != fX[Nm1] && fX[Nm1] < *(mBinBoundaries[Nm1].begin() + indexBounds.second))
0188 tmpIt -= 1;
0189 fN[Nm1] = *tmpIt;
0190 }
0191
0192 tuple_type to_find =
0193 gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY>([&](size_t i) { return (i < SIZE) ? fN[i] : -9999; });
0194 return (mIndexMap.find(to_find) != mIndexMap.end()) ? mIndexMap.at(to_find) : -1;
0195 }