Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-16 03:19:40

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // JetResolution
0004 // -------------
0005 //
0006 //            11/05/2010 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
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 // GLOBAL FUNCTION DEFINITIONS
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 // CONSTRUCTION / DESTRUCTION
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 // IMPLEMENTATION OF MEMBER FUNCTIONS
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 // IMPLEMENTATION OF GLOBAL FUNCTIONS
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 }