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