JetMETObjects/src/SimpleJetCorrectionUncertainty.cc

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)  // true = UP
0048       value.push_back(p[ind + 1]);
0049     else  // false = DOWN
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   // Linear interpolation through the points (x[i],y[i]). First find the line that
0070   // is defined by the points and then calculate the y(z).
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 }