File indexing completed on 2025-01-04 00:29:26
0001 #include "CondFormats/JetMETObjects/interface/SimpleJetCorrectionUncertainty.h"
0002 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include <vector>
0005 #include <string>
0006
0007
0008 SimpleJetCorrectionUncertainty::SimpleJetCorrectionUncertainty(const std::string& fDataFile) : mParameters{fDataFile} {}
0009
0010 SimpleJetCorrectionUncertainty::SimpleJetCorrectionUncertainty(const JetCorrectorParameters& fParameters)
0011 : mParameters{fParameters} {}
0012
0013 float SimpleJetCorrectionUncertainty::uncertainty(const std::vector<float>& fX, float fY, bool fDirection) const {
0014 float result = 1.;
0015 int bin = mParameters.binIndex(fX);
0016 if (bin < 0) {
0017 edm::LogError("SimpleJetCorrectionUncertainty") << " bin variables out of range";
0018 result = -999.0;
0019 } else
0020 result = uncertaintyBin((unsigned)bin, fY, fDirection);
0021 return result;
0022 }
0023
0024 namespace {
0025 class Span3 {
0026 public:
0027 Span3(std::vector<float> const& iFrom, int offset = 0) : start_{iFrom.data() + offset}, size_{iFrom.size() / 3} {}
0028 float last() const { return fromIndex(size_ - 1); }
0029 float operator[](std::size_t i) const { return fromIndex(i); }
0030 std::size_t size() const { return size_; }
0031
0032 private:
0033 float fromIndex(std::size_t i) const { return *(start_ + 3 * i); }
0034 float const* start_;
0035 std::size_t size_;
0036 };
0037 int findBin(const Span3& v, float x) {
0038 int n = v.size() - 1;
0039 if (n <= 0)
0040 return -1;
0041 if (x < v[0] || x >= v[n])
0042 return -1;
0043 for (int i = 0; i < n; i++) {
0044 if (x >= v[i] && x < v[i + 1])
0045 return i;
0046 }
0047 return 0;
0048 }
0049 }
0050 float SimpleJetCorrectionUncertainty::uncertaintyBin(unsigned fBin, float fY, bool fDirection) const {
0051 if (fBin >= mParameters.size()) {
0052 edm::LogError("SimpleJetCorrectionUncertainty")
0053 << " wrong bin: " << fBin << ": only " << mParameters.size() << " are available";
0054 return -999.0;
0055 }
0056 const std::vector<float>& p = mParameters.record(fBin).parameters();
0057 if ((p.size() % 3) != 0 or p.empty())
0058 throw cms::Exception("SimpleJetCorrectionUncertainty")
0059 << "wrong # of parameters: multiple of 3 expected, " << p.size() << " got";
0060 Span3 yGrid{p};
0061 Span3 value{p, fDirection ? 1 : 2};
0062 float result = -1.0;
0063 if (fY <= yGrid[0])
0064 result = value[0];
0065 else if (fY >= yGrid.last())
0066 result = value.last();
0067 else {
0068 int bin = findBin(yGrid, fY);
0069 float vx[2], vy[2];
0070 for (int i = 0; i < 2; i++) {
0071 vx[i] = yGrid[bin + i];
0072 vy[i] = value[bin + i];
0073 }
0074 result = linearInterpolation(fY, vx, vy);
0075 }
0076 return result;
0077 }
0078
0079 float SimpleJetCorrectionUncertainty::linearInterpolation(float fZ, const float fX[2], const float fY[2]) const {
0080
0081
0082 float r = 0;
0083 if (fX[0] == fX[1]) {
0084 if (fY[0] == fY[1])
0085 r = fY[0];
0086 else {
0087 edm::LogError("SimpleJetCorrectionUncertainty") << " interpolation error";
0088 return -999.0;
0089 }
0090 } else {
0091 float a = (fY[1] - fY[0]) / (fX[1] - fX[0]);
0092 float b = (fY[0] * fX[1] - fY[1] * fX[0]) / (fX[1] - fX[0]);
0093 r = a * fZ + b;
0094 }
0095 return r;
0096 }