JetMETObjects/src/JetCorrectorParametersHelper.cc

0001 //
0002 // Original Author:  Alexx Perloff Feb 22, 2017
0003 // $Id: JetCorrectorParameters.cc,v 1.20 2012/03/01 18:24:53 srappocc Exp $
0004 //
0005 // Helper class for JetCorrectorParameters
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 //--- JetCorrectorParameters::JetCorrectorParametersHelper initializer ---
0016 //--- initializes the mBinMap for quick lookup of mRecords ---------------
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           //mMap.emplace(gen_tuple<SIZE-1>([&](size_t k){return mRecordsLocal[i].xMin(k);}),std::make_pair(start,end));
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   // This function checks that all of the middle binned parameters (not the first or last) contain the same range of bins.
0085   // If not an exception will be thrown as the mapping will work only if all bins but the last are uniform.
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   //The order of looping (records or dimensions) does not matter because you have to go through all of them once anyway
0096   //Loop over each binned dimension that isn't the first or the last
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     //Loop over the mRecords vector
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 //--- JetCorrectorParameters::JetCorrectorParametersHelper sanity checks -
0122 //--- checks that some conditions are met before finding the bin index ---
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 //--- JetCorrectorParameters::JetCorrectorParametersHelper binIndexN -----
0142 //--- returns the index of the record defined by fX (non-linear search) --
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   //Create a container for the indices
0150   std::vector<float> fN(SIZE, -1);
0151   std::vector<float>::const_iterator tmpIt;
0152
0153   // make sure that fX are within the first and last boundaries of mBinBoundaries (other than last dimension)
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     // lower_bound finds the entry with the next highest value to fX[0]
0159     // so unless the two values are equal, you want the next lowest bin boundary
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   //find the index bounds for the possible values of the last dimension
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     //Check that the requested value is within the bin boundaries for the last dimension
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 }