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(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 }  // namespace
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   // Linear interpolation through the points (x[i],y[i]). First find the line that
0081   // is defined by the points and then calculate the y(z).
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 }