00001
00002
00003
00004
00005
00006
00008
00009
00010 #include "CondFormats/JetMETObjects/interface/JetResolution.h"
00011 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
00012
00013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00014
00015 #include <TMath.h>
00016
00017
00018 #include <iostream>
00019 #include <sstream>
00020 #include <cassert>
00021
00022
00023 using namespace std;
00024
00025
00027
00029
00030
00031 double fnc_dscb(double*xx,double*pp);
00032 double fnc_gaussalpha(double*xx,double*pp);
00033 double fnc_gaussalpha1alpha2(double*xx,double*pp);
00034
00035
00037
00039
00040
00041 JetResolution::JetResolution()
00042 : resolutionFnc_(0)
00043 {
00044 resolutionFnc_ = new TF1();
00045 }
00046
00047
00048
00049 JetResolution::JetResolution(const string& fileName,bool doGaussian)
00050 : resolutionFnc_(0)
00051 {
00052 initialize(fileName,doGaussian);
00053 }
00054
00055
00056
00057 JetResolution::~JetResolution()
00058 {
00059 delete resolutionFnc_;
00060 for (unsigned i=0;i<parameterFncs_.size();i++) delete parameterFncs_[i];
00061 for (unsigned i=0;i<parameters_.size();i++) delete parameters_[i];
00062 }
00063
00064
00066
00068
00069
00070 void JetResolution::initialize(const string& fileName,bool doGaussian)
00071 {
00072 size_t pos;
00073
00074 name_ = fileName;
00075 pos = name_.find_last_of('.'); name_ = name_.substr(0,pos);
00076 pos = name_.find_last_of('/'); name_ = name_.substr(pos+1);
00077
00078 JetCorrectorParameters resolutionPars(fileName,"resolution");
00079 string fncname = "fResolution_" + name_;
00080 string formula = resolutionPars.definitions().formula();
00081 if (doGaussian) resolutionFnc_=new TF1(fncname.c_str(),"gaus",0.,5.);
00082 else if (formula=="DSCB") resolutionFnc_=new TF1(fncname.c_str(),fnc_dscb,0.,5.,7);
00083 else if (formula=="GaussAlpha1Alpha2") resolutionFnc_=new TF1(fncname.c_str(),fnc_gaussalpha1alpha2,-5.,5.,5);
00084 else if (formula=="GaussAlpha") resolutionFnc_=new TF1(fncname.c_str(),fnc_gaussalpha,-5.,5.,4);
00085 else resolutionFnc_=new TF1(fncname.c_str(),formula.c_str(),0.,5.);
00086
00087 resolutionFnc_->SetNpx(200);
00088 resolutionFnc_->SetParName(0,"N");
00089 resolutionFnc_->SetParameter(0,1.0);
00090 unsigned nPar(1);
00091
00092 string tmp = resolutionPars.definitions().level();
00093 pos = tmp.find(':');
00094 while (!tmp.empty()) {
00095 string paramAsStr = tmp.substr(0,pos);
00096 if (!doGaussian||paramAsStr=="mean"||paramAsStr=="sigma") {
00097 parameters_.push_back(new JetCorrectorParameters(fileName,paramAsStr));
00098 formula = parameters_.back()->definitions().formula();
00099 parameterFncs_.push_back(new TF1(("f"+paramAsStr+"_"+name()).c_str(),formula.c_str(),
00100 parameters_.back()->record(0).parameters()[0],
00101 parameters_.back()->record(0).parameters()[1]));
00102 resolutionFnc_->SetParName(nPar,parameters_.back()->definitions().level().c_str());
00103 nPar++;
00104 }
00105 tmp = (pos==string::npos) ? "" : tmp.substr(pos+1);
00106 pos = tmp.find(':');
00107 }
00108
00109 assert(nPar==(unsigned)resolutionFnc_->GetNpar());
00110 assert(!doGaussian||nPar==3);
00111 }
00112
00113
00114
00115 TF1* JetResolution::resolutionEtaPt(float eta, float pt) const
00116 {
00117 vector<float> x; x.push_back(eta);
00118 vector<float> y; y.push_back(pt);
00119 return resolution(x,y);
00120 }
00121
00122
00123
00124 TF1* JetResolution::resolution(const vector<float>& x,
00125 const vector<float>& y) const
00126 {
00127 unsigned N(y.size());
00128 for (unsigned iPar=0;iPar<parameters_.size();iPar++) {
00129 int bin = parameters_[iPar]->binIndex(x);
00130 assert(bin>=0);
00131 assert(bin<(int)parameters_[iPar]->size());
00132 const std::vector<float>& pars = parameters_[iPar]->record(bin).parameters();
00133 for (unsigned i=2*N;i<pars.size();i++)
00134 parameterFncs_[iPar]->SetParameter(i-2*N,pars[i]);
00135 float yy[4];
00136 for (unsigned i=0;i<N;i++)
00137 yy[i] = (y[i] < pars[2*i]) ? pars[2*i] : (y[i] > pars[2*i+1]) ? pars[2*i+1] : y[i];
00138 resolutionFnc_->SetParameter(iPar+1,
00139 parameterFncs_[iPar]->Eval(yy[0],yy[1],yy[2],yy[3]));
00140 }
00141 return resolutionFnc_;
00142 }
00143
00144
00145
00146 TF1* JetResolution::parameterEta(const string& parameterName, float eta)
00147 {
00148 vector<float> x; x.push_back(eta);
00149 return parameter(parameterName,x);
00150 }
00151
00152
00153
00154 TF1* JetResolution::parameter(const string& parameterName,const vector<float>& x)
00155 {
00156 TF1* result(0);
00157 for (unsigned i=0;i<parameterFncs_.size()&&result==0;i++) {
00158 string fncname = parameterFncs_[i]->GetName();
00159 if (fncname.find("f"+parameterName)==0) {
00160 stringstream ssname; ssname<<parameterFncs_[i]->GetName();
00161 for (unsigned ii=0;ii<x.size();ii++)
00162 ssname<<"_"<<parameters_[i]->definitions().binVar(ii)<<x[ii];
00163 result = (TF1*)parameterFncs_[i]->Clone();
00164 result->SetName(ssname.str().c_str());
00165 int N = parameters_[i]->definitions().nParVar();
00166 int bin = parameters_[i]->binIndex(x);
00167 assert(bin>=0);
00168 assert(bin<(int)parameters_[i]->size());
00169 const std::vector<float>& pars = parameters_[i]->record(bin).parameters();
00170 for (unsigned ii=2*N;ii<pars.size();ii++) result->SetParameter(ii-2*N,pars[ii]);
00171 }
00172 }
00173
00174 if (0==result) cerr<<"JetResolution::parameter() ERROR: no parameter "
00175 <<parameterName<<" found."<<endl;
00176
00177 return result;
00178 }
00179
00180
00181
00182 double JetResolution::parameterEtaEval(const std::string& parameterName, float eta, float pt)
00183 {
00184 TF1* func(0);
00185 JetCorrectorParameters* params(0);
00186 for (std::vector<TF1*>::size_type ifunc = 0; ifunc < parameterFncs_.size(); ++ifunc)
00187 {
00188 std::string fncname = parameterFncs_[ifunc]->GetName();
00189 if ( !(fncname.find("f"+parameterName) == 0) ) continue;
00190 params = parameters_[ifunc];
00191 func = (TF1*)parameterFncs_[ifunc];
00192 break;
00193 }
00194
00195 if (!func)
00196 edm::LogError("ParameterNotFound") << "JetResolution::parameterEtaEval(): no parameter \""
00197 << parameterName << "\" found" << std::endl;
00198
00199 std::vector<float> etas; etas.push_back(eta);
00200 int bin = params->binIndex(etas);
00201
00202 if ( !(0 <= bin && bin < (int)params->size() ) )
00203 edm::LogError("ParameterNotFound") << "JetResolution::parameterEtaEval(): bin out of range: "
00204 << bin << std::endl;
00205
00206 const std::vector<float>& pars = params->record(bin).parameters();
00207
00208 int N = params->definitions().nParVar();
00209 for (unsigned ii = 2*N; ii < pars.size(); ++ii)
00210 {
00211 func->SetParameter(ii-2*N, pars[ii]);
00212 }
00213
00214 return func->Eval(pt);
00215 }
00216
00217
00219
00221
00222
00223 double fnc_dscb(double*xx,double*pp)
00224 {
00225 double x = xx[0];
00226 double N = pp[0];
00227 double mu = pp[1];
00228 double sig = pp[2];
00229 double a1 = pp[3];
00230 double p1 = pp[4];
00231 double a2 = pp[5];
00232 double p2 = pp[6];
00233
00234 double u = (x-mu)/sig;
00235 double A1 = TMath::Power(p1/TMath::Abs(a1),p1)*TMath::Exp(-a1*a1/2);
00236 double A2 = TMath::Power(p2/TMath::Abs(a2),p2)*TMath::Exp(-a2*a2/2);
00237 double B1 = p1/TMath::Abs(a1) - TMath::Abs(a1);
00238 double B2 = p2/TMath::Abs(a2) - TMath::Abs(a2);
00239
00240 double result(N);
00241 if (u<-a1) result *= A1*TMath::Power(B1-u,-p1);
00242 else if (u<a2) result *= TMath::Exp(-u*u/2);
00243 else result *= A2*TMath::Power(B2+u,-p2);
00244 return result;
00245 }
00246
00247
00248
00249 double fnc_gaussalpha(double *v, double *par)
00250 {
00251 double N =par[0];
00252 double mean =par[1];
00253 double sigma=par[2];
00254 double alpha=par[3];
00255 double t =TMath::Abs((v[0]-mean)/sigma);
00256 double cut = 1.0;
00257 return (t<=cut) ? N*TMath::Exp(-0.5*t*t) : N*TMath::Exp(-0.5*(alpha*(t-cut)+cut*cut));
00258 }
00259
00260
00261
00262 double fnc_gaussalpha1alpha2(double *v, double *par)
00263 {
00264 double N =par[0];
00265 double mean =par[1];
00266 double sigma =par[2];
00267 double alpha1=par[3];
00268 double alpha2=par[4];
00269 double t =TMath::Abs((v[0]-mean)/sigma);
00270 double cut = 1.0;
00271 return
00272 (t<=cut) ? N*TMath::Exp(-0.5*t*t) :
00273 ((v[0]-mean)>=0) ? N*TMath::Exp(-0.5*(alpha1*(t-cut)+cut*cut)) :
00274 N*TMath::Exp(-0.5*(alpha2*(t-cut)+cut*cut));
00275 }
00276
1.7.3