1 #ifndef PhysicsTools_Heppy_MuScleFitCorrector_Functions_h
2 #define PhysicsTools_Heppy_MuScleFitCorrector_Functions_h
26 ParSet(
const TString & inputName,
const double & inputStep,
const double & inputMini,
const double & inputMaxi) :
31 std::cout <<
"setting name = " << inputName << std::endl;
46 virtual double scale(
const double &
pt,
const double &
eta,
const double &
phi,
const int chg,
const T & parScale)
const = 0;
51 virtual void setPar(
double* Start,
double*
Step,
double* Mini,
double* Maxi,
int* ind,
52 TString* parname,
const T & parResol,
const std::vector<int> & parResolOrder,
const std::vector<ParSet> & parSet ) {
53 if(
int(parSet.size()) != this->
parNum_ ) {
54 std::cout <<
"Error: wrong number of parameter initializations = " << parSet.size() <<
". Number of parameters is " << this->
parNum_ << std::endl;
57 for(
int iPar=0; iPar<this->
parNum_; ++iPar ) {
58 Start[iPar] = parResol[iPar];
59 Step[iPar] = parSet[iPar].step;
60 Mini[iPar] = parSet[iPar].mini;
61 Maxi[iPar] = parSet[iPar].maxi;
62 ind[iPar] = parResolOrder[iPar];
63 parname[iPar] = parSet[iPar].name;
79 virtual double scale(
const double &
pt,
const double &
eta,
const double &
phi,
const int chg,
const T & parScale)
const {
80 double ampl(0), phase(0), twist(0), ampl2(0), freq2(0), phase2(0);
83 if ( eta < parScale[4] ) {
84 ampl = parScale[1]; phase = parScale[2]; ampl2 = parScale[21]; freq2 = parScale[22]; phase2 = parScale[23];
85 twist = parScale[3]*(eta-parScale[4])+parScale[7]*(parScale[4]-parScale[8])+parScale[11]*parScale[8];
87 }
else if ( parScale[4] <= eta && eta < parScale[8] ) {
88 ampl = parScale[5]; phase = parScale[6];
89 twist = parScale[7]*(eta-parScale[8])+parScale[11]*parScale[8] ;
91 }
else if ( parScale[8] <= eta && eta < parScale[12] ) {
92 ampl = parScale[9]; phase = parScale[10];
93 twist = parScale[11]*
eta;
95 }
else if ( parScale[12] <= eta && eta < parScale[16] ) {
96 ampl = parScale[13]; phase = parScale[14];
97 twist = parScale[15]*(eta-parScale[12])+parScale[11]*parScale[12];
99 }
else if ( parScale[16] < eta ) {
100 ampl = parScale[17]; phase = parScale[18]; ampl2 = parScale[24]; freq2 = parScale[25]; phase2 = parScale[26];
101 twist = parScale[19]*(eta-parScale[16])+parScale[15]*(parScale[16]-parScale[12])+parScale[11]*parScale[12];
105 double curv = (1.+parScale[0])*((
double)chg/pt
108 -ampl2*
sin(freq2*phi+phase2)
110 return 1./((double)chg*curv);
123 virtual double sigmaPt(
const double &
pt,
const double &
eta,
const T & parval) = 0;
140 inline double getGEO(
const double &
pt,
const double &
eta,
const T & parval){
144 inline double getMS(
const double &
pt,
const double &
eta,
const T & parval){
145 if( eta < -2.0 )
return( parval[1] );
146 if( eta < -1.8 )
return( parval[2] );
147 if( eta < -1.6 )
return( parval[3] );
148 if( eta < -1.2 )
return( parval[4] );
149 if( eta < -0.8 )
return( parval[5] );
150 if( eta < 0. )
return( parval[6] );
151 if( eta < 0.8 )
return( parval[7] );
152 if( eta < 1.2 )
return( parval[8] );
153 if( eta < 1.6 )
return( parval[9] );
154 if( eta < 1.8 )
return( parval[10] );
155 if( eta < 2.0 )
return( parval[11] );
156 return( parval[12] );
159 virtual double sigmaPt(
const double &
pt,
const double &
eta,
const T & parval)
161 return pt*
getGEO(pt,eta,parval) +
getMS(pt,eta,parval);
176 if( eta < -2.0 )
return 1;
177 if( eta < -1.8 )
return 2;
178 if( eta < -1.6 )
return 3;
179 if( eta < -1.2 )
return 4;
180 if( eta < -0.8 )
return 5;
181 if( eta < 0. )
return 6;
182 if( eta < 0.8 )
return 7;
183 if( eta < 1.2 )
return 8;
184 if( eta < 1.6 )
return 9;
185 if( eta < 1.8 )
return 10;
186 if( eta < 2.0 )
return 11;
190 virtual double sigmaPt(
const double &
pt,
const double &
eta,
const T & parval)
205 inline double getGEO(
const double &
pt,
const double &
eta,
const T & parval){
208 if ( eta < parval[0] ) {
209 qGEO = parval[10]*(eta-parval[0])*(eta-parval[0])+parval[11];
210 }
else if ( parval[0] <= eta && eta < parval[3] ) {
213 qGEO = parval[12]*(eta-parval[3])*(eta-parval[3])+parval[11];
219 return parval[4] + parval[5]*eta*
eta;
223 return parval[15] + parval[16]*eta*
eta;
227 return parval[6] + parval[7]*(eta-parval[0])*(eta-parval[0]);
231 return parval[8] + parval[9]*(eta-parval[3])*(eta-parval[3]);
235 double x_1 = parval[13];
237 double x_2 = parval[0];
239 return( (eta - x_1)*(y_2 - y_1)/(x_2 - x_1) + y_1 );
243 double x_1 = parval[14];
245 double x_2 = parval[3];
247 return( (eta - x_1)*(y_2 - y_1)/(x_2 - x_1) + y_1 );
251 inline double getMSC(
const double &
pt,
const double &
eta,
const T & parval){
254 if ( eta < parval[0] ) {
256 }
else if ( parval[0] <= eta && eta < parval[13] ) {
258 }
else if ( parval[13] <= eta && eta < parval[1] ) {
260 }
else if ( parval[1] <= eta && eta < parval[2] ) {
262 }
else if ( parval[2] <= eta && eta < parval[14] ) {
264 }
else if ( parval[14] <= eta && eta < parval[3] ) {
272 virtual double sigmaPt(
const double &
pt,
const double &
eta,
const T & parval)
274 double qGEO =
getGEO(pt, eta, parval);
275 double qMSC =
getMSC(pt, eta, parval);
284 switch ( identifier ) {
286 default:
std::cout <<
"scaleFunctionService error: wrong identifier = " << identifier << std::endl;
exit(1);
293 switch ( identifier ) {
297 default:
std::cout <<
"resolutionFunctService error: wrong identifier = " << identifier << std::endl;
exit(1);
scaleFunctionBase< double * > * scaleFunctionService(const int identifier)
Service to build the scale functor corresponding to the passed identifier.
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)
double centralParabola(const double &pt, const double &eta, const T &parval)
double rightParabola(const double &pt, const double &eta, const T &parval)
double getMSC(const double &pt, const double &eta, const T &parval)
virtual double scale(const double &pt, const double &eta, const double &phi, const int chg, const T &parScale) const
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)
Sin< T >::type sin(const T &t)
double leftLine(const double &pt, const double &eta, const T &parval)
int etaBin(const double &eta)
virtual ~scaleFunctionBase()=0
double middleParabola(const double &pt, const double &eta, const T &parval)
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)=0
virtual double sigmaPt(const double &pt, const double &eta, const T &parval)
double getMS(const double &pt, const double &eta, const T &parval)
virtual int parNum() const
ParSet(const TString &inputName, const double &inputStep, const double &inputMini, const double &inputMaxi)
virtual double scale(const double &pt, const double &eta, const double &phi, const int chg, const T &parScale) const =0
double getGEO(const double &pt, const double &eta, const T &parval)
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
virtual void setPar(double *Start, double *Step, double *Mini, double *Maxi, int *ind, TString *parname, const T &parResol, const std::vector< int > &parResolOrder, const std::vector< ParSet > &parSet)
double rightLine(const double &pt, const double &eta, const T &parval)
virtual int parNum() const
virtual ~resolutionFunctionBase()=0
Power< A, B >::type pow(const A &a, const B &b)
double getGEO(const double &pt, const double &eta, const T &parval)
double leftParabola(const double &pt, const double &eta, const T &parval)