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