JetMETObjects/src/SimpleJetCorrector.cc

0001 #include "CondFormats/JetMETObjects/interface/SimpleJetCorrector.h"
0002 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0003 #include "CondFormats/JetMETObjects/interface/Utilities.h"
0004 #include <iostream>
0005 #include <sstream>
0006 #include <cmath>
0007 #include <algorithm>
0008
0009 //------------------------------------------------------------------------
0010 //--- SimpleJetCorrector constructor -------------------------------------
0011 //--- reads arguments from a file ----------------------------------------
0012 //------------------------------------------------------------------------
0013 SimpleJetCorrector::SimpleJetCorrector(const std::string& fDataFile, const std::string& fOption)
0014     : mParameters(fDataFile, fOption), mFunc((mParameters.definitions()).formula()) {
0015   mDoInterpolation = false;
0016   if (mParameters.definitions().isResponse())
0017     mInvertVar = findInvertVar();
0018 }
0019 //------------------------------------------------------------------------
0020 //--- SimpleJetCorrector constructor -------------------------------------
0021 //--- reads arguments from a file ----------------------------------------
0022 //------------------------------------------------------------------------
0023 SimpleJetCorrector::SimpleJetCorrector(const JetCorrectorParameters& fParameters)
0024     : mParameters(fParameters), mFunc((mParameters.definitions()).formula()) {
0025   mDoInterpolation = false;
0026   if (mParameters.definitions().isResponse())
0027     mInvertVar = findInvertVar();
0028 }
0029
0030 //------------------------------------------------------------------------
0031 //--- calculates the correction ------------------------------------------
0032 //------------------------------------------------------------------------
0033 float SimpleJetCorrector::correction(const std::vector<float>& fX, const std::vector<float>& fY) const {
0034   float result = 1.;
0035   float tmp = 0.0;
0036   float cor = 0.0;
0037   int bin = -1;
0038   bin = (fX.size() <= 3 && !fX.empty()) ? mParameters.binIndexN(fX) : mParameters.binIndex(fX);
0039   if (bin < 0)
0040     return result;
0041   if (!mDoInterpolation)
0042     result = correctionBin(bin, fY);
0043   else {
0044     for (unsigned i = 0; i < mParameters.definitions().nBinVar(); i++) {
0045       float xMiddle[3];
0046       float xValue[3];
0047       int prevBin = mParameters.neighbourBin((unsigned)bin, i, false);
0048       int nextBin = mParameters.neighbourBin((unsigned)bin, i, true);
0049       if (prevBin >= 0 && nextBin >= 0) {
0050         xMiddle[0] = mParameters.record(prevBin).xMiddle(i);
0051         xMiddle[1] = mParameters.record(bin).xMiddle(i);
0052         xMiddle[2] = mParameters.record(nextBin).xMiddle(i);
0053         xValue[0] = correctionBin(prevBin, fY);
0054         xValue[1] = correctionBin(bin, fY);
0055         xValue[2] = correctionBin(nextBin, fY);
0056         cor = quadraticInterpolation(fX[i], xMiddle, xValue);
0057         tmp += cor;
0058       } else {
0059         cor = correctionBin(bin, fY);
0060         tmp += cor;
0061       }
0062     }
0063     result = tmp / mParameters.definitions().nBinVar();
0064   }
0065   if (result <= 0) {
0066     edm::LogWarning("SimpleJetCorrector")
0067         << "Null or negative jet energy correction factor evaluated: " << result << ". Truncating to 10e-10.";
0068     result = 10e-10;
0069   }
0070   return result;
0071 }
0072 //------------------------------------------------------------------------
0073 //--- calculates the correction for a specific bin -----------------------
0074 //------------------------------------------------------------------------
0075 float SimpleJetCorrector::correctionBin(unsigned fBin, const std::vector<float>& fY) const {
0076   if (fBin >= mParameters.size()) {
0077     std::stringstream sserr;
0078     sserr << "wrong bin: " << fBin << ": only " << mParameters.size() << " available!";
0079     handleError("SimpleJetCorrector", sserr.str());
0080   }
0081   unsigned N = fY.size();
0082   if (!fY.empty() && N <= 4) {
0083     const std::vector<float>& par = mParameters.record(fBin).parameters();
0084     auto nParams = static_cast<int>(par.size() - 2 * N);
0085     if (nParams > 0) {
0086       double params[nParams];
0087       for (unsigned int i = 2 * N; i < par.size(); i++) {
0088         params[i - 2 * N] = par[i];
0089       }
0090       double x[4] = {};
0091       for (unsigned i = 0; i < N; i++) {
0092         x[i] = (fY[i] < par[2 * i]) ? par[2 * i] : (fY[i] > par[2 * i + 1]) ? par[2 * i + 1] : fY[i];
0093       }
0094       if (mParameters.definitions().isResponse()) {
0095         return invert(x, params);
0096       }
0097       return mFunc.evaluate(reco::formula::ArrayAdaptor(x, N), reco::formula::ArrayAdaptor(params, par.size() - 2 * N));
0098     } else {
0099       return 1.0;
0100     }
0101   } else {
0102     edm::LogWarning("SimpleJetCorrector") << "Expect number of variables to be 0<N<=4 but is " << N << "in bin " << fBin
0103                                           << ". Returning zero" << std::endl;
0104     return 0.0;
0105   }
0106 }
0107 //------------------------------------------------------------------------
0108 //--- find invertion variable (JetPt) ------------------------------------
0109 //------------------------------------------------------------------------
0110 unsigned SimpleJetCorrector::findInvertVar() {
0111   unsigned result = 9999;
0112   std::vector<std::string> vv = mParameters.definitions().parVar();
0113   for (unsigned i = 0; i < vv.size(); i++)
0114     if (vv[i] == "JetPt") {
0115       result = i;
0116       break;
0117     }
0118   if (result >= vv.size())
0119     handleError("SimpleJetCorrector", "Response inversion is required but JetPt is not specified as parameter");
0120   return result;
0121 }
0122 //------------------------------------------------------------------------
0123 //--- inversion ----------------------------------------------------------
0124 //------------------------------------------------------------------------
0125 float SimpleJetCorrector::invert(const double* args, const double* params) const {
0126   unsigned nMax = 50;
0127   float precision = 0.0001;
0128   float rsp = 1.0;
0129   float e = 1.0;
0130   double x[4];
0131   unsigned nLoop = 0;
0132
0133   // 4 dimensions (x, y, z, t) supported in TFormula
0134   memcpy(&x, args, sizeof(double) * 4);
0135
0136   while (e > precision && nLoop < nMax) {
0137     rsp = mFunc.evaluate(reco::formula::ArrayAdaptor(x, 4),
0138                          reco::formula::ArrayAdaptor(params, mFunc.numberOfParameters()));
0139     float tmp = x[mInvertVar] * rsp;
0140     e = fabs(tmp - args[mInvertVar]) / args[mInvertVar];
0141     x[mInvertVar] = args[mInvertVar] / rsp;
0142     nLoop++;
0143   }
0144   return 1. / rsp;
0145 }