JetMETObjects/src/JetResolution.cc

////////////////////////////////////////////////////////////////////////////////
//
// JetResolution
// -------------
//
//            11/05/2010 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
////////////////////////////////////////////////////////////////////////////////


#include "CondFormats/JetMETObjects/interface/JetResolution.h"
#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"


#include <TMath.h>


#include <iostream>
#include <sstream>
#include <cassert>


using namespace std;


////////////////////////////////////////////////////////////////////////////////
// GLOBAL FUNCTION DEFINITIONS
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
double fnc_dscb(double*xx,double*pp);
double fnc_gaussalpha(double*xx,double*pp);
double fnc_gaussalpha1alpha2(double*xx,double*pp);


////////////////////////////////////////////////////////////////////////////////
// CONSTRUCTION / DESTRUCTION
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
JetResolution::JetResolution()
  : resolutionFnc_(0)
{
  resolutionFnc_ = new TF1();
}


//______________________________________________________________________________
JetResolution::JetResolution(const string& fileName,bool doGaussian)
  : resolutionFnc_(0)
{
  initialize(fileName,doGaussian);
}


//______________________________________________________________________________
JetResolution::~JetResolution()
{
  delete resolutionFnc_;
  for (unsigned i=0;i<parameterFncs_.size();i++) delete parameterFncs_[i];
  for (unsigned i=0;i<parameters_.size();i++)    delete parameters_[i];
}


////////////////////////////////////////////////////////////////////////////////
// IMPLEMENTATION OF MEMBER FUNCTIONS
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
void JetResolution::initialize(const string& fileName,bool doGaussian)
{
  size_t pos;
  
  name_ = fileName;
  pos = name_.find_last_of('.'); name_ = name_.substr(0,pos);
  pos = name_.find_last_of('/'); name_ = name_.substr(pos+1);
  
  JetCorrectorParameters resolutionPars(fileName,"resolution");
  string fncname = "fResolution_" + name_;
  string formula = resolutionPars.definitions().formula();
  if      (doGaussian)                   resolutionFnc_=new TF1(fncname.c_str(),"gaus",0.,5.);
  else if (formula=="DSCB")              resolutionFnc_=new TF1(fncname.c_str(),fnc_dscb,0.,5.,7);
  else if (formula=="GaussAlpha1Alpha2") resolutionFnc_=new TF1(fncname.c_str(),fnc_gaussalpha1alpha2,-5.,5.,5);
  else if (formula=="GaussAlpha")        resolutionFnc_=new TF1(fncname.c_str(),fnc_gaussalpha,-5.,5.,4);
  else                                   resolutionFnc_=new TF1(fncname.c_str(),formula.c_str(),0.,5.);
  
  resolutionFnc_->SetNpx(200);
  resolutionFnc_->SetParName(0,"N");
  resolutionFnc_->SetParameter(0,1.0);
  unsigned nPar(1);
  
  string tmp = resolutionPars.definitions().level();
  pos = tmp.find(':');
  while (!tmp.empty()) {
    string paramAsStr = tmp.substr(0,pos);
    if (!doGaussian||paramAsStr=="mean"||paramAsStr=="sigma") {
      parameters_.push_back(new JetCorrectorParameters(fileName,paramAsStr));
      formula = parameters_.back()->definitions().formula();
      parameterFncs_.push_back(new TF1(("f"+paramAsStr+"_"+name()).c_str(),formula.c_str(),
                                       parameters_.back()->record(0).parameters()[0],
                                       parameters_.back()->record(0).parameters()[1]));
      resolutionFnc_->SetParName(nPar,parameters_.back()->definitions().level().c_str());
      nPar++;
    }
    tmp = (pos==string::npos) ? "" : tmp.substr(pos+1);
    pos = tmp.find(':');
  }
  
  assert(nPar==(unsigned)resolutionFnc_->GetNpar());
  assert(!doGaussian||nPar==3);
}
  

//______________________________________________________________________________
TF1* JetResolution::resolutionEtaPt(float eta, float pt) const
{
  vector<float> x; x.push_back(eta);
  vector<float> y; y.push_back(pt);
  return resolution(x,y);
}


//______________________________________________________________________________
TF1* JetResolution::resolution(const vector<float>& x,
                               const vector<float>& y) const
{
  unsigned N(y.size());
  for (unsigned iPar=0;iPar<parameters_.size();iPar++) {
    int bin = parameters_[iPar]->binIndex(x);
    assert(bin>=0);
    assert(bin<(int)parameters_[iPar]->size());
    const std::vector<float>& pars = parameters_[iPar]->record(bin).parameters();
    for (unsigned i=2*N;i<pars.size();i++)
      parameterFncs_[iPar]->SetParameter(i-2*N,pars[i]);
    float yy[4];
    for (unsigned i=0;i<N;i++)
      yy[i] = (y[i] < pars[2*i]) ? pars[2*i] : (y[i] > pars[2*i+1]) ? pars[2*i+1] : y[i];
    resolutionFnc_->SetParameter(iPar+1,
                                 parameterFncs_[iPar]->Eval(yy[0],yy[1],yy[2],yy[3]));
  }
  return resolutionFnc_;
}


//______________________________________________________________________________
TF1* JetResolution::parameterEta(const string& parameterName, float eta)
{
  vector<float> x; x.push_back(eta);
  return parameter(parameterName,x);
}


//______________________________________________________________________________
TF1* JetResolution::parameter(const string& parameterName,const vector<float>& x)
{
  TF1* result(0);
  for (unsigned i=0;i<parameterFncs_.size()&&result==0;i++) {
    string fncname = parameterFncs_[i]->GetName();
    if (fncname.find("f"+parameterName)==0) {
      stringstream ssname; ssname<<parameterFncs_[i]->GetName();
      for (unsigned ii=0;ii<x.size();ii++)
        ssname<<"_"<<parameters_[i]->definitions().binVar(ii)<<x[ii];
      result = (TF1*)parameterFncs_[i]->Clone();
      result->SetName(ssname.str().c_str());
      int N = parameters_[i]->definitions().nParVar();
      int bin = parameters_[i]->binIndex(x);
      assert(bin>=0);
      assert(bin<(int)parameters_[i]->size());
      const std::vector<float>& pars = parameters_[i]->record(bin).parameters();
      for (unsigned ii=2*N;ii<pars.size();ii++) result->SetParameter(ii-2*N,pars[ii]); 
    }
  }
  
  if (0==result) cerr<<"JetResolution::parameter() ERROR: no parameter "
                     <<parameterName<<" found."<<endl;
  
  return result;
}


////////////////////////////////////////////////////////////////////////////////
// IMPLEMENTATION OF GLOBAL FUNCTIONS
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
double fnc_dscb(double*xx,double*pp)
{
  double x   = xx[0];
  double N   = pp[0];
  double mu  = pp[1];
  double sig = pp[2];
  double a1  = pp[3];
  double p1  = pp[4];
  double a2  = pp[5];
  double p2  = pp[6];
  
  double u   = (x-mu)/sig;
  double A1  = TMath::Power(p1/TMath::Abs(a1),p1)*TMath::Exp(-a1*a1/2);
  double A2  = TMath::Power(p2/TMath::Abs(a2),p2)*TMath::Exp(-a2*a2/2);
  double B1  = p1/TMath::Abs(a1) - TMath::Abs(a1);
  double B2  = p2/TMath::Abs(a2) - TMath::Abs(a2);

  double result(N);
  if      (u<-a1) result *= A1*TMath::Power(B1-u,-p1);
  else if (u<a2)  result *= TMath::Exp(-u*u/2);
  else            result *= A2*TMath::Power(B2+u,-p2);
  return result;
}


//______________________________________________________________________________
double fnc_gaussalpha(double *v, double *par)
{
    double N    =par[0];
    double mean =par[1];
    double sigma=par[2];
    double alpha=par[3];
    double t    =TMath::Abs((v[0]-mean)/sigma);
    double cut  = 1.0;
    return (t<=cut) ? N*TMath::Exp(-0.5*t*t) : N*TMath::Exp(-0.5*(alpha*(t-cut)+cut*cut));
}


//______________________________________________________________________________
double fnc_gaussalpha1alpha2(double *v, double *par)
{
    double N     =par[0];
    double mean  =par[1];
    double sigma =par[2];
    double alpha1=par[3];
    double alpha2=par[4];
    double t     =TMath::Abs((v[0]-mean)/sigma);
    double cut = 1.0;
    return
      (t<=cut) ? N*TMath::Exp(-0.5*t*t) :
      ((v[0]-mean)>=0) ? N*TMath::Exp(-0.5*(alpha1*(t-cut)+cut*cut)) :
      N*TMath::Exp(-0.5*(alpha2*(t-cut)+cut*cut));
}

Revision:	1.1
Committed:	Wed May 30 09:20:33 2012 UTC (12 years, 11 months ago) by peiffer
Content type:	text/plain
Branch:	MAIN
CVS Tags:	v1-00, Feb-15-2013-v1, Feb-14-2013, Feb-07-2013-v1, Jan-17-2013-v2, Jan-17-2013-v1, Jan-16-2012-v1, Jan-09-2012-v2, Jan-09-2012-v1, Dec-26-2012-v1, Dec-20-2012-v1, Dec-17-2012-v1, Nov-30-2012-v2, Nov-30-2012-v1, HEAD
Log Message:	added JetMETObjects
#	User	Rev	Content
1	peiffer	1.1	////////////////////////////////////////////////////////////////////////////////
2			//
3			// JetResolution
4			// -------------
5			//
6			// 11/05/2010 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
7			////////////////////////////////////////////////////////////////////////////////
8
9
10			#include "CondFormats/JetMETObjects/interface/JetResolution.h"
11			#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
12
13
14			#include <TMath.h>
15
16
17			#include <iostream>
18			#include <sstream>
19			#include <cassert>
20
21
22			using namespace std;
23
24
25			////////////////////////////////////////////////////////////////////////////////
26			// GLOBAL FUNCTION DEFINITIONS
27			////////////////////////////////////////////////////////////////////////////////
28
29			//______________________________________________________________________________
30			double fnc_dscb(doublexx,doublepp);
31			double fnc_gaussalpha(doublexx,doublepp);
32			double fnc_gaussalpha1alpha2(doublexx,doublepp);
33
34
35			////////////////////////////////////////////////////////////////////////////////
36			// CONSTRUCTION / DESTRUCTION
37			////////////////////////////////////////////////////////////////////////////////
38
39			//______________________________________________________________________________
40			JetResolution::JetResolution()
41			: resolutionFnc_(0)
42			{
43			resolutionFnc_ = new TF1();
44			}
45
46
47			//______________________________________________________________________________
48			JetResolution::JetResolution(const string& fileName,bool doGaussian)
49			: resolutionFnc_(0)
50			{
51			initialize(fileName,doGaussian);
52			}
53
54
55			//______________________________________________________________________________
56			JetResolution::~JetResolution()
57			{
58			delete resolutionFnc_;
59			for (unsigned i=0;i<parameterFncs_.size();i++) delete parameterFncs_[i];
60			for (unsigned i=0;i<parameters_.size();i++) delete parameters_[i];
61			}
62
63
64			////////////////////////////////////////////////////////////////////////////////
65			// IMPLEMENTATION OF MEMBER FUNCTIONS
66			////////////////////////////////////////////////////////////////////////////////
67
68			//______________________________________________________________________________
69			void JetResolution::initialize(const string& fileName,bool doGaussian)
70			{
71			size_t pos;
72
73			name_ = fileName;
74			pos = name_.find_last_of('.'); name_ = name_.substr(0,pos);
75			pos = name_.find_last_of('/'); name_ = name_.substr(pos+1);
76
77			JetCorrectorParameters resolutionPars(fileName,"resolution");
78			string fncname = "fResolution_" + name_;
79			string formula = resolutionPars.definitions().formula();
80			if (doGaussian) resolutionFnc_=new TF1(fncname.c_str(),"gaus",0.,5.);
81			else if (formula=="DSCB") resolutionFnc_=new TF1(fncname.c_str(),fnc_dscb,0.,5.,7);
82			else if (formula=="GaussAlpha1Alpha2") resolutionFnc_=new TF1(fncname.c_str(),fnc_gaussalpha1alpha2,-5.,5.,5);
83			else if (formula=="GaussAlpha") resolutionFnc_=new TF1(fncname.c_str(),fnc_gaussalpha,-5.,5.,4);
84			else resolutionFnc_=new TF1(fncname.c_str(),formula.c_str(),0.,5.);
85
86			resolutionFnc_->SetNpx(200);
87			resolutionFnc_->SetParName(0,"N");
88			resolutionFnc_->SetParameter(0,1.0);
89			unsigned nPar(1);
90
91			string tmp = resolutionPars.definitions().level();
92			pos = tmp.find(':');
93			while (!tmp.empty()) {
94			string paramAsStr = tmp.substr(0,pos);
95			if (!doGaussian\|\|paramAsStr=="mean"\|\|paramAsStr=="sigma") {
96			parameters_.push_back(new JetCorrectorParameters(fileName,paramAsStr));
97			formula = parameters_.back()->definitions().formula();
98			parameterFncs_.push_back(new TF1(("f"+paramAsStr+"_"+name()).c_str(),formula.c_str(),
99			parameters_.back()->record(0).parameters()[0],
100			parameters_.back()->record(0).parameters()[1]));
101			resolutionFnc_->SetParName(nPar,parameters_.back()->definitions().level().c_str());
102			nPar++;
103			}
104			tmp = (pos==string::npos) ? "" : tmp.substr(pos+1);
105			pos = tmp.find(':');
106			}
107
108			assert(nPar==(unsigned)resolutionFnc_->GetNpar());
109			assert(!doGaussian\|\|nPar==3);
110			}
111
112
113			//______________________________________________________________________________
114			TF1* JetResolution::resolutionEtaPt(float eta, float pt) const
115			{
116			vector<float> x; x.push_back(eta);
117			vector<float> y; y.push_back(pt);
118			return resolution(x,y);
119			}
120
121
122			//______________________________________________________________________________
123			TF1* JetResolution::resolution(const vector<float>& x,
124			const vector<float>& y) const
125			{
126			unsigned N(y.size());
127			for (unsigned iPar=0;iPar<parameters_.size();iPar++) {
128			int bin = parameters_[iPar]->binIndex(x);
129			assert(bin>=0);
130			assert(bin<(int)parameters_[iPar]->size());
131			const std::vector<float>& pars = parameters_[iPar]->record(bin).parameters();
132			for (unsigned i=2*N;i<pars.size();i++)
133			parameterFncs_[iPar]->SetParameter(i-2*N,pars[i]);
134			float yy[4];
135			for (unsigned i=0;i<N;i++)
136			yy[i] = (y[i] < pars[2i]) ? pars[2i] : (y[i] > pars[2i+1]) ? pars[2i+1] : y[i];
137			resolutionFnc_->SetParameter(iPar+1,
138			parameterFncs_[iPar]->Eval(yy[0],yy[1],yy[2],yy[3]));
139			}
140			return resolutionFnc_;
141			}
142
143
144			//______________________________________________________________________________
145			TF1* JetResolution::parameterEta(const string& parameterName, float eta)
146			{
147			vector<float> x; x.push_back(eta);
148			return parameter(parameterName,x);
149			}
150
151
152			//______________________________________________________________________________
153			TF1* JetResolution::parameter(const string& parameterName,const vector<float>& x)
154			{
155			TF1* result(0);
156			for (unsigned i=0;i<parameterFncs_.size()&&result==0;i++) {
157			string fncname = parameterFncs_[i]->GetName();
158			if (fncname.find("f"+parameterName)==0) {
159			stringstream ssname; ssname<<parameterFncs_[i]->GetName();
160			for (unsigned ii=0;ii<x.size();ii++)
161			ssname<<"_"<<parameters_[i]->definitions().binVar(ii)<<x[ii];
162			result = (TF1*)parameterFncs_[i]->Clone();
163			result->SetName(ssname.str().c_str());
164			int N = parameters_[i]->definitions().nParVar();
165			int bin = parameters_[i]->binIndex(x);
166			assert(bin>=0);
167			assert(bin<(int)parameters_[i]->size());
168			const std::vector<float>& pars = parameters_[i]->record(bin).parameters();
169			for (unsigned ii=2N;ii<pars.size();ii++) result->SetParameter(ii-2N,pars[ii]);
170			}
171			}
172
173			if (0==result) cerr<<"JetResolution::parameter() ERROR: no parameter "
174			<<parameterName<<" found."<<endl;
175
176			return result;
177			}
178
179
180			////////////////////////////////////////////////////////////////////////////////
181			// IMPLEMENTATION OF GLOBAL FUNCTIONS
182			////////////////////////////////////////////////////////////////////////////////
183
184			//______________________________________________________________________________
185			double fnc_dscb(doublexx,doublepp)
186			{
187			double x = xx[0];
188			double N = pp[0];
189			double mu = pp[1];
190			double sig = pp[2];
191			double a1 = pp[3];
192			double p1 = pp[4];
193			double a2 = pp[5];
194			double p2 = pp[6];
195
196			double u = (x-mu)/sig;
197			double A1 = TMath::Power(p1/TMath::Abs(a1),p1)TMath::Exp(-a1a1/2);
198			double A2 = TMath::Power(p2/TMath::Abs(a2),p2)TMath::Exp(-a2a2/2);
199			double B1 = p1/TMath::Abs(a1) - TMath::Abs(a1);
200			double B2 = p2/TMath::Abs(a2) - TMath::Abs(a2);
201
202			double result(N);
203			if (u<-a1) result = A1TMath::Power(B1-u,-p1);
204			else if (u<a2) result = TMath::Exp(-uu/2);
205			else result = A2TMath::Power(B2+u,-p2);
206			return result;
207			}
208
209
210			//______________________________________________________________________________
211			double fnc_gaussalpha(double v, double par)
212			{
213			double N =par[0];
214			double mean =par[1];
215			double sigma=par[2];
216			double alpha=par[3];
217			double t =TMath::Abs((v[0]-mean)/sigma);
218			double cut = 1.0;
219			return (t<=cut) ? NTMath::Exp(-0.5tt) : NTMath::Exp(-0.5(alpha(t-cut)+cut*cut));
220			}
221
222
223			//______________________________________________________________________________
224			double fnc_gaussalpha1alpha2(double v, double par)
225			{
226			double N =par[0];
227			double mean =par[1];
228			double sigma =par[2];
229			double alpha1=par[3];
230			double alpha2=par[4];
231			double t =TMath::Abs((v[0]-mean)/sigma);
232			double cut = 1.0;
233			return
234			(t<=cut) ? NTMath::Exp(-0.5t*t) :
235			((v[0]-mean)>=0) ? NTMath::Exp(-0.5(alpha1(t-cut)+cutcut)) :
236			NTMath::Exp(-0.5(alpha2(t-cut)+cutcut));
237			}
238