File indexing completed on 2023-03-17 10:47:08
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() { mParameters = new JetCorrectorParameters(); }
0009
0010 SimpleJetCorrectionUncertainty::SimpleJetCorrectionUncertainty(const std::string& fDataFile) {
0011 mParameters = new JetCorrectorParameters(fDataFile);
0012 }
0013
0014 SimpleJetCorrectionUncertainty::SimpleJetCorrectionUncertainty(const JetCorrectorParameters& fParameters) {
0015 mParameters = new JetCorrectorParameters(fParameters);
0016 }
0017
0018 SimpleJetCorrectionUncertainty::~SimpleJetCorrectionUncertainty() { delete mParameters; }
0019
0020 float SimpleJetCorrectionUncertainty::uncertainty(const std::vector<float>& fX, float fY, bool fDirection) const {
0021 float result = 1.;
0022 int bin = mParameters->binIndex(fX);
0023 if (bin < 0) {
0024 edm::LogError("SimpleJetCorrectionUncertainty") << " bin variables out of range";
0025 result = -999.0;
0026 } else
0027 result = uncertaintyBin((unsigned)bin, fY, fDirection);
0028 return result;
0029 }
0030
0031 float SimpleJetCorrectionUncertainty::uncertaintyBin(unsigned fBin, float fY, bool fDirection) const {
0032 if (fBin >= mParameters->size()) {
0033 edm::LogError("SimpleJetCorrectionUncertainty")
0034 << " wrong bin: " << fBin << ": only " << mParameters->size() << " are available";
0035 return -999.0;
0036 }
0037 const std::vector<float>& p = mParameters->record(fBin).parameters();
0038 if ((p.size() % 3) != 0)
0039 throw cms::Exception("SimpleJetCorrectionUncertainty")
0040 << "wrong # of parameters: multiple of 3 expected, " << p.size() << " got";
0041 std::vector<float> yGrid, value;
0042 unsigned int N = p.size() / 3;
0043 float result = -1.0;
0044 for (unsigned i = 0; i < N; i++) {
0045 unsigned ind = 3 * i;
0046 yGrid.push_back(p[ind]);
0047 if (fDirection)
0048 value.push_back(p[ind + 1]);
0049 else
0050 value.push_back(p[ind + 2]);
0051 }
0052 if (fY <= yGrid[0])
0053 result = value[0];
0054 else if (fY >= yGrid[N - 1])
0055 result = value[N - 1];
0056 else {
0057 int bin = findBin(yGrid, fY);
0058 float vx[2], vy[2];
0059 for (int i = 0; i < 2; i++) {
0060 vx[i] = yGrid[bin + i];
0061 vy[i] = value[bin + i];
0062 }
0063 result = linearInterpolation(fY, vx, vy);
0064 }
0065 return result;
0066 }
0067
0068 float SimpleJetCorrectionUncertainty::linearInterpolation(float fZ, const float fX[2], const float fY[2]) const {
0069
0070
0071 float r = 0;
0072 if (fX[0] == fX[1]) {
0073 if (fY[0] == fY[1])
0074 r = fY[0];
0075 else {
0076 edm::LogError("SimpleJetCorrectionUncertainty") << " interpolation error";
0077 return -999.0;
0078 }
0079 } else {
0080 float a = (fY[1] - fY[0]) / (fX[1] - fX[0]);
0081 float b = (fY[0] * fX[1] - fY[1] * fX[0]) / (fX[1] - fX[0]);
0082 r = a * fZ + b;
0083 }
0084 return r;
0085 }
0086
0087 int SimpleJetCorrectionUncertainty::findBin(const std::vector<float>& v, float x) const {
0088 int i;
0089 int n = v.size() - 1;
0090 if (n <= 0)
0091 return -1;
0092 if (x < v[0] || x >= v[n])
0093 return -1;
0094 for (i = 0; i < n; i++) {
0095 if (x >= v[i] && x < v[i + 1])
0096 return i;
0097 }
0098 return 0;
0099 }