File indexing completed on 2021-02-14 12:51:16
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
0010
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
0020
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
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
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
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
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
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 }