File indexing completed on 2024-04-06 12:02:18
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "CondFormats/JetMETObjects/interface/JetResolution.h"
0010 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include <TMath.h>
0015
0016 #include <iostream>
0017 #include <sstream>
0018 #include <cassert>
0019
0020 using namespace std;
0021
0022
0023
0024
0025
0026
0027 double fnc_dscb(double* xx, double* pp);
0028 double fnc_gaussalpha(double* xx, double* pp);
0029 double fnc_gaussalpha1alpha2(double* xx, double* pp);
0030
0031
0032
0033
0034
0035
0036 JetResolution::JetResolution() : resolutionFnc_(nullptr) { resolutionFnc_ = new TF1(); }
0037
0038
0039 JetResolution::JetResolution(const string& fileName, bool doGaussian) : resolutionFnc_(nullptr) {
0040 initialize(fileName, doGaussian);
0041 }
0042
0043
0044 JetResolution::~JetResolution() {
0045 delete resolutionFnc_;
0046 for (unsigned i = 0; i < parameterFncs_.size(); i++)
0047 delete parameterFncs_[i];
0048 for (unsigned i = 0; i < parameters_.size(); i++)
0049 delete parameters_[i];
0050 }
0051
0052
0053
0054
0055
0056
0057 void JetResolution::initialize(const string& fileName, bool doGaussian) {
0058 size_t pos;
0059
0060 name_ = fileName;
0061 pos = name_.find_last_of('.');
0062 name_ = name_.substr(0, pos);
0063 pos = name_.find_last_of('/');
0064 name_ = name_.substr(pos + 1);
0065
0066 JetCorrectorParameters resolutionPars(fileName, "resolution");
0067 string fncname = "fResolution_" + name_;
0068 string formula = resolutionPars.definitions().formula();
0069 if (doGaussian)
0070 resolutionFnc_ = new TF1(fncname.c_str(), "gaus", 0., 5.);
0071 else if (formula == "DSCB")
0072 resolutionFnc_ = new TF1(fncname.c_str(), fnc_dscb, 0., 5., 7);
0073 else if (formula == "GaussAlpha1Alpha2")
0074 resolutionFnc_ = new TF1(fncname.c_str(), fnc_gaussalpha1alpha2, -5., 5., 5);
0075 else if (formula == "GaussAlpha")
0076 resolutionFnc_ = new TF1(fncname.c_str(), fnc_gaussalpha, -5., 5., 4);
0077 else
0078 resolutionFnc_ = new TF1(fncname.c_str(), formula.c_str(), 0., 5.);
0079
0080 resolutionFnc_->SetNpx(200);
0081 resolutionFnc_->SetParName(0, "N");
0082 resolutionFnc_->SetParameter(0, 1.0);
0083 unsigned nPar(1);
0084
0085 string tmp = resolutionPars.definitions().level();
0086 pos = tmp.find(':');
0087 while (!tmp.empty()) {
0088 string paramAsStr = tmp.substr(0, pos);
0089 if (!doGaussian || paramAsStr == "mean" || paramAsStr == "sigma") {
0090 parameters_.push_back(new JetCorrectorParameters(fileName, paramAsStr));
0091 formula = parameters_.back()->definitions().formula();
0092 parameterFncs_.push_back(new TF1(("f" + paramAsStr + "_" + name()).c_str(),
0093 formula.c_str(),
0094 parameters_.back()->record(0).parameters()[0],
0095 parameters_.back()->record(0).parameters()[1]));
0096 resolutionFnc_->SetParName(nPar, parameters_.back()->definitions().level().c_str());
0097 nPar++;
0098 }
0099 tmp = (pos == string::npos) ? "" : tmp.substr(pos + 1);
0100 pos = tmp.find(':');
0101 }
0102
0103 assert(nPar == (unsigned)resolutionFnc_->GetNpar());
0104 assert(!doGaussian || nPar == 3);
0105 }
0106
0107
0108 TF1* JetResolution::resolutionEtaPt(float eta, float pt) const {
0109 vector<float> x;
0110 x.push_back(eta);
0111 vector<float> y;
0112 y.push_back(pt);
0113 return resolution(x, y);
0114 }
0115
0116
0117 TF1* JetResolution::resolution(const vector<float>& x, const vector<float>& y) const {
0118 unsigned N(y.size());
0119 for (unsigned iPar = 0; iPar < parameters_.size(); iPar++) {
0120 int bin = parameters_[iPar]->binIndex(x);
0121 assert(bin >= 0);
0122 assert(bin < (int)parameters_[iPar]->size());
0123 const std::vector<float>& pars = parameters_[iPar]->record(bin).parameters();
0124 for (unsigned i = 2 * N; i < pars.size(); i++)
0125 parameterFncs_[iPar]->SetParameter(i - 2 * N, pars[i]);
0126 float yy[4] = {};
0127 for (unsigned i = 0; i < N; i++)
0128 yy[i] = (y[i] < pars[2 * i]) ? pars[2 * i] : (y[i] > pars[2 * i + 1]) ? pars[2 * i + 1] : y[i];
0129 resolutionFnc_->SetParameter(iPar + 1, parameterFncs_[iPar]->Eval(yy[0], yy[1], yy[2], yy[3]));
0130 }
0131 return resolutionFnc_;
0132 }
0133
0134
0135 TF1* JetResolution::parameterEta(const string& parameterName, float eta) {
0136 vector<float> x;
0137 x.push_back(eta);
0138 return parameter(parameterName, x);
0139 }
0140
0141
0142 TF1* JetResolution::parameter(const string& parameterName, const vector<float>& x) {
0143 TF1* result(nullptr);
0144 for (unsigned i = 0; i < parameterFncs_.size() && result == nullptr; i++) {
0145 string fncname = parameterFncs_[i]->GetName();
0146 if (fncname.find("f" + parameterName) == 0) {
0147 stringstream ssname;
0148 ssname << parameterFncs_[i]->GetName();
0149 for (unsigned ii = 0; ii < x.size(); ii++)
0150 ssname << "_" << parameters_[i]->definitions().binVar(ii) << x[ii];
0151 result = (TF1*)parameterFncs_[i]->Clone();
0152 result->SetName(ssname.str().c_str());
0153 int N = parameters_[i]->definitions().nParVar();
0154 int bin = parameters_[i]->binIndex(x);
0155 assert(bin >= 0);
0156 assert(bin < (int)parameters_[i]->size());
0157 const std::vector<float>& pars = parameters_[i]->record(bin).parameters();
0158 for (unsigned ii = 2 * N; ii < pars.size(); ii++)
0159 result->SetParameter(ii - 2 * N, pars[ii]);
0160 }
0161 }
0162
0163 if (nullptr == result)
0164 cerr << "JetResolution::parameter() ERROR: no parameter " << parameterName << " found." << endl;
0165
0166 return result;
0167 }
0168
0169
0170 double JetResolution::parameterEtaEval(const std::string& parameterName, float eta, float pt) {
0171 TF1* func(nullptr);
0172 JetCorrectorParameters* params(nullptr);
0173 for (std::vector<TF1*>::size_type ifunc = 0; ifunc < parameterFncs_.size(); ++ifunc) {
0174 std::string fncname = parameterFncs_[ifunc]->GetName();
0175 if (!(fncname.find("f" + parameterName) == 0))
0176 continue;
0177 params = parameters_[ifunc];
0178 func = (TF1*)parameterFncs_[ifunc];
0179 break;
0180 }
0181
0182 if (!func)
0183 edm::LogError("ParameterNotFound") << "JetResolution::parameterEtaEval(): no parameter \"" << parameterName
0184 << "\" found" << std::endl;
0185
0186 std::vector<float> etas;
0187 etas.push_back(eta);
0188 int bin = params->binIndex(etas);
0189
0190 if (!(0 <= bin && bin < (int)params->size()))
0191 edm::LogError("ParameterNotFound") << "JetResolution::parameterEtaEval(): bin out of range: " << bin << std::endl;
0192
0193 const std::vector<float>& pars = params->record(bin).parameters();
0194
0195 int N = params->definitions().nParVar();
0196 for (unsigned ii = 2 * N; ii < pars.size(); ++ii) {
0197 func->SetParameter(ii - 2 * N, pars[ii]);
0198 }
0199
0200 return func->Eval(pt);
0201 }
0202
0203
0204
0205
0206
0207
0208 double fnc_dscb(double* xx, double* pp) {
0209 double x = xx[0];
0210 double N = pp[0];
0211 double mu = pp[1];
0212 double sig = pp[2];
0213 double a1 = pp[3];
0214 double p1 = pp[4];
0215 double a2 = pp[5];
0216 double p2 = pp[6];
0217
0218 double u = (x - mu) / sig;
0219 double A1 = TMath::Power(p1 / TMath::Abs(a1), p1) * TMath::Exp(-a1 * a1 / 2);
0220 double A2 = TMath::Power(p2 / TMath::Abs(a2), p2) * TMath::Exp(-a2 * a2 / 2);
0221 double B1 = p1 / TMath::Abs(a1) - TMath::Abs(a1);
0222 double B2 = p2 / TMath::Abs(a2) - TMath::Abs(a2);
0223
0224 double result(N);
0225 if (u < -a1)
0226 result *= A1 * TMath::Power(B1 - u, -p1);
0227 else if (u < a2)
0228 result *= TMath::Exp(-u * u / 2);
0229 else
0230 result *= A2 * TMath::Power(B2 + u, -p2);
0231 return result;
0232 }
0233
0234
0235 double fnc_gaussalpha(double* v, double* par) {
0236 double N = par[0];
0237 double mean = par[1];
0238 double sigma = par[2];
0239 double alpha = par[3];
0240 double t = TMath::Abs((v[0] - mean) / sigma);
0241 double cut = 1.0;
0242 return (t <= cut) ? N * TMath::Exp(-0.5 * t * t) : N * TMath::Exp(-0.5 * (alpha * (t - cut) + cut * cut));
0243 }
0244
0245
0246 double fnc_gaussalpha1alpha2(double* v, double* par) {
0247 double N = par[0];
0248 double mean = par[1];
0249 double sigma = par[2];
0250 double alpha1 = par[3];
0251 double alpha2 = par[4];
0252 double t = TMath::Abs((v[0] - mean) / sigma);
0253 double cut = 1.0;
0254 return (t <= cut) ? N * TMath::Exp(-0.5 * t * t)
0255 : ((v[0] - mean) >= 0) ? N * TMath::Exp(-0.5 * (alpha1 * (t - cut) + cut * cut))
0256 : N * TMath::Exp(-0.5 * (alpha2 * (t - cut) + cut * cut));
0257 }
