00001 #ifndef FUNCTIONS_H
00002 #define FUNCTIONS_H
00003
00004 #include <iostream>
00005 #include <vector>
00006 #include <cmath>
00007 #include "TMath.h"
00008 #include "TString.h"
00009 #include "TF1.h"
00010 #include "TRandom.h"
00011 #include "MuonAnalysis/MomentumScaleCalibration/interface/SigmaPtDiff.h"
00012
00016 struct ParameterSet
00017 {
00018 ParameterSet() {}
00019 ParameterSet(const TString & inputName, const double & inputStep, const double & inputMini, const double & inputMaxi) :
00020 step(inputStep),
00021 mini(inputMini),
00022 maxi(inputMaxi)
00023 {
00024 std::cout << "setting name = " << inputName << std::endl;
00025 name = inputName;
00026 }
00027 TString name;
00028 double step, mini, maxi;
00029 };
00030
00031
00032
00033
00043 template <class T>
00044 class scaleFunctionBase {
00045 public:
00046 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const = 0;
00047 virtual ~scaleFunctionBase() = 0;
00049 virtual void resetParameters(std::vector<double> * scaleVec) const {
00050 std::cout << "ERROR: the resetParameters method must be defined in each scale function" << std::endl;
00051 std::cout << "Please add it to the scaleFunction you are using" << std::endl;
00052 exit(1);
00053 }
00055 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) = 0;
00056 virtual int parNum() const { return parNum_; }
00057 protected:
00058 int parNum_;
00060 virtual void setPar(double* Start, double* Step, double* Mini, double* Maxi, int* ind,
00061 TString* parname, const T & parScale, const std::vector<int> & parScaleOrder,
00062 double* thisStep, double* thisMini, double* thisMaxi, TString* thisParName ) {
00063 for( int iPar=0; iPar<this->parNum_; ++iPar ) {
00064 Start[iPar] = parScale[iPar];
00065 Step[iPar] = thisStep[iPar];
00066 Mini[iPar] = thisMini[iPar];
00067 Maxi[iPar] = thisMaxi[iPar];
00068 ind[iPar] = parScaleOrder[iPar];
00069 parname[iPar] = thisParName[iPar];
00070 }
00071 }
00072 };
00073 template <class T> inline scaleFunctionBase<T>::~scaleFunctionBase() { }
00074
00075
00076 template <class T>
00077 class scaleFunctionType0 : public scaleFunctionBase<T> {
00078 public:
00079 scaleFunctionType0() {
00080
00081
00082 this->parNum_ = 0;
00083 }
00084 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const { return pt; }
00085 virtual void resetParameters(std::vector<double> * scaleVec) const {}
00086 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {}
00087 };
00088
00089
00090 template <class T>
00091 class scaleFunctionType1 : public scaleFunctionBase<T> {
00092 public:
00093 scaleFunctionType1() { this->parNum_ = 2; }
00094 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00095 return ( (parScale[0] + parScale[1]*pt)*pt );
00096 }
00097
00098 virtual void resetParameters(std::vector<double> * scaleVec) const {
00099 scaleVec->push_back(1);
00100 for( int i=1; i<this->parNum_; ++i ) {
00101 scaleVec->push_back(0);
00102 }
00103 }
00104 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00105 double thisStep[] = {0.001, 0.01};
00106 TString thisParName[] = {"Pt offset", "Pt slope"};
00107 if( muonType == 1 ) {
00108 double thisMini[] = {0.97, -0.1};
00109 double thisMaxi[] = {1.03, 0.1};
00110 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00111 } else {
00112 double thisMini[] = {0.97, -0.1};
00113 double thisMaxi[] = {1.03, 0.1};
00114 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00115 }
00116 }
00117 };
00118
00119
00120 template <class T>
00121 class scaleFunctionType2 : public scaleFunctionBase<T> {
00122 public:
00123 scaleFunctionType2() { this->parNum_ = 2; }
00124 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00125 return ( (parScale[0] + parScale[1]*std::fabs(eta))*pt );
00126 }
00127
00128 virtual void resetParameters(std::vector<double> * scaleVec) const {
00129 scaleVec->push_back(1);
00130 for( int i=1; i<this->parNum_; ++i ) {
00131 scaleVec->push_back(0);
00132 }
00133 }
00134 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00135 double thisStep[] = {0.001, 0.01};
00136 TString thisParName[] = {"Eta offset", "Eta slope"};
00137 if( muonType == 1 ) {
00138 double thisMini[] = {0.9, -0.3};
00139 double thisMaxi[] = {1.1, 0.3};
00140 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00141 } else {
00142 double thisMini[] = {0.97, -0.1};
00143 double thisMaxi[] = {1.03, 0.1};
00144 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00145 }
00146 }
00147 };
00148
00149
00150 template <class T>
00151 class scaleFunctionType3 : public scaleFunctionBase<T> {
00152 public:
00153 scaleFunctionType3() { this->parNum_ = 4; }
00154 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00155 return( (parScale[0] + parScale[1]*sin(parScale[2]*phi + parScale[3]))*pt );
00156 }
00157
00158 virtual void resetParameters(std::vector<double> * scaleVec) const {
00159 scaleVec->push_back(1);
00160 for( int i=1; i<this->parNum_; ++i ) {
00161 scaleVec->push_back(0);
00162 }
00163 }
00164 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00165 double thisStep[] = {0.001, 0.01, 0.01, 0.01};
00166 TString thisParName[] = {"Phi offset", "Phi ampl", "Phi freq", "Phi phase"};
00167 if( muonType == 1 ) {
00168 double thisMini[] = {0.9, -0.1, -0.1, -0.1};
00169 double thisMaxi[] = {1.1, 0.1, 0.1, 0.1};
00170 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00171 } else {
00172 double thisMini[] = {0.97, -0.05, 6, -3.14};
00173 double thisMaxi[] = {1.03, 0.05, 0.1, 3.14};
00174 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00175 }
00176 }
00177 };
00178
00179
00180 template <class T>
00181 class scaleFunctionType4 : public scaleFunctionBase<T> {
00182 public:
00183 scaleFunctionType4() { this->parNum_ = 3; }
00184 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00185 return( (parScale[0] + parScale[1]*pt +
00186 parScale[2]*std::fabs(eta))*pt );
00187 }
00188
00189 virtual void resetParameters(std::vector<double> * scaleVec) const {
00190 scaleVec->push_back(1);
00191 for( int i=1; i<this->parNum_; ++i ) {
00192 scaleVec->push_back(0);
00193 }
00194 }
00195 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00196 double thisStep[] = {0.001, 0.01, 0.01};
00197 TString thisParName[] = {"Pt offset", "Pt slope", "Eta slope"};
00198 if( muonType == 1 ) {
00199 double thisMini[] = {0.9, -0.1, -0.1};
00200 double thisMaxi[] = {1.1, 0.1, 0.1};
00201 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00202 } else {
00203 double thisMini[] = {0.97, -0.02, -0.02};
00204 double thisMaxi[] = {1.03, 0.02, 0.02};
00205 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00206 }
00207 }
00208 };
00209
00210
00211 template <class T>
00212 class scaleFunctionType5 : public scaleFunctionBase<T> {
00213 public:
00214 scaleFunctionType5() { this->parNum_ = 3; }
00215 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00216 return( (parScale[0] + parScale[1]*pt +
00217 parScale[2]*sin(phi))*pt );
00218 }
00219
00220 virtual void resetParameters(std::vector<double> * scaleVec) const {
00221 scaleVec->push_back(1);
00222 for( int i=1; i<this->parNum_; ++i ) {
00223 scaleVec->push_back(0);
00224 }
00225 }
00226 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00227 double thisStep[] = {0.001, 0.01, 0.01};
00228 TString thisParName[] = {"Pt offset", "Pt slope", "Phi ampl"};
00229 if( muonType == 1 ) {
00230 double thisMini[] = {0.9, -0.1, -0.3};
00231 double thisMaxi[] = {1.1, 0.1, 0.3};
00232 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00233 } else {
00234 double thisMini[] = {0.97, -0.02, -0.3};
00235 double thisMaxi[] = {1.03, 0.02, 0.3};
00236 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00237 }
00238 }
00239 };
00240
00241
00242 template <class T>
00243 class scaleFunctionType6 : public scaleFunctionBase<T> {
00244 public:
00245 scaleFunctionType6() { this->parNum_ = 3; }
00246 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00247 return( (parScale[0] + parScale[1]*std::fabs(eta) +
00248 parScale[2]*sin(phi))*pt );
00249 }
00250
00251 virtual void resetParameters(std::vector<double> * scaleVec) const {
00252 scaleVec->push_back(1);
00253 for( int i=1; i<this->parNum_; ++i ) {
00254 scaleVec->push_back(0);
00255 }
00256 }
00257 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind,
00258 TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00259 double thisStep[] = {0.001, 0.01, 0.01};
00260 TString thisParName[] = {"Eta offset", "Eta slope", "Phi ampl"};
00261 if( muonType == 1 ) {
00262 double thisMini[] = {0.9, -0.1, -0.3};
00263 double thisMaxi[] = {1.1, 0.1, 0.3};
00264 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00265 } else {
00266 double thisMini[] = {0.97, -0.02, -0.3};
00267 double thisMaxi[] = {1.03, 0.02, 0.3};
00268 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00269 }
00270 }
00271 };
00272
00273
00274 template <class T>
00275 class scaleFunctionType7 : public scaleFunctionBase<T> {
00276 public:
00277 scaleFunctionType7() { this->parNum_ = 4; }
00278 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00279 return( (parScale[0] + parScale[1]*pt +
00280 parScale[2]*std::fabs(eta) +
00281 parScale[3]*sin(phi))*pt );
00282 }
00283
00284 virtual void resetParameters(std::vector<double> * scaleVec) const {
00285 scaleVec->push_back(1);
00286 for( int i=1; i<this->parNum_; ++i ) {
00287 scaleVec->push_back(0);
00288 }
00289 }
00290 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind,
00291 TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00292 double thisStep[] = {0.001, 0.01, 0.01, 0.01};
00293 TString thisParName[] = {"Pt offset", "Pt slope", "Eta slope", "Phi ampl"};
00294 if( muonType == 1 ) {
00295 double thisMini[] = {0.9, -0.1, -0.3, -0.3};
00296 double thisMaxi[] = {1.1, 0.1, 0.3, 0.3};
00297 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00298 } else {
00299 double thisMini[] = {0.97, -0.02, -0.3, -0.3};
00300 double thisMaxi[] = {1.03, 0.02, 0.3, 0.3};
00301 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00302 }
00303 }
00304 };
00305
00306
00307 template <class T>
00308 class scaleFunctionType8 : public scaleFunctionBase<T> {
00309 public:
00310 scaleFunctionType8() { this->parNum_ = 4; }
00311 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00312 return( (parScale[0] + parScale[1]*pt +
00313 parScale[2]*std::fabs(eta) +
00314 parScale[3]*eta*eta)*pt );
00315 }
00316
00317 virtual void resetParameters(std::vector<double> * scaleVec) const {
00318 scaleVec->push_back(1);
00319 for( int i=1; i<this->parNum_; ++i ) {
00320 scaleVec->push_back(0);
00321 }
00322 }
00323 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00324 double thisStep[] = {0.00001, 0.000001, 0.0000001, 0.0000001};
00325 TString thisParName[] = {"Pt offset", "Pt slope", "Eta slope", "Eta quadr"};
00326 if( muonType == 1 ) {
00327 double thisMini[] = {0.9, -0.3, -0.3, -0.3};
00328 double thisMaxi[] = {1.1, 0.3, 0.3, 0.3};
00329 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00330 } else {
00331
00332
00333 double thisMini[] = {0.9, -0.002, -0.01, -0.005};
00334 double thisMaxi[] = {1.1, 0.002, 0.01, 0.005};
00335 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00336 }
00337 }
00338 };
00339
00340
00341 template <class T>
00342 class scaleFunctionType9 : public scaleFunctionBase<T> {
00343 public:
00344 scaleFunctionType9() { this->parNum_ = 2; }
00345 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00346 return( (parScale[0] + exp(parScale[1]*pt))*pt );
00347 }
00348
00349 virtual void resetParameters(std::vector<double> * scaleVec) const {
00350 scaleVec->push_back(1);
00351 for( int i=1; i<this->parNum_; ++i ) {
00352 scaleVec->push_back(0);
00353 }
00354 }
00355 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00356 double thisStep[] = {0.001, 0.01};
00357 TString thisParName[] = {"Pt offset", "Pt expon"};
00358 double thisMini[] = {0.97, -0.1};
00359 double thisMaxi[] = {1.03, 0.1};
00360 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00361 }
00362 };
00363
00364
00365 template <class T>
00366 class scaleFunctionType10 : public scaleFunctionBase<T> {
00367 public:
00368 scaleFunctionType10() { this->parNum_ = 3; }
00369 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00370 return( (parScale[0] + parScale[1]*pt +
00371 parScale[2]*pt*pt)*pt );
00372 }
00373
00374 virtual void resetParameters(std::vector<double> * scaleVec) const {
00375 scaleVec->push_back(1);
00376 for( int i=1; i<this->parNum_; ++i ) {
00377 scaleVec->push_back(0);
00378 }
00379 }
00380 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00381 double thisStep[] = {0.001, 0.01, 0.01};
00382 TString thisParName[] = {"Pt offset", "Pt slope", "Pt quadr"};
00383 double thisMini[] = {0.97, -0.1, -0.001};
00384 double thisMaxi[] = {1.03, 0.1, 0.001};
00385 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00386 }
00387 };
00388
00389
00390 template <class T>
00391 class scaleFunctionType11 : public scaleFunctionBase<T> {
00392 public:
00393 scaleFunctionType11() { this->parNum_ = 4; }
00394 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00395 return( (parScale[0] + parScale[1]*pt +
00396 (double)chg*parScale[2]*sin(phi+parScale[3]))*pt );
00397 }
00398
00399 virtual void resetParameters(std::vector<double> * scaleVec) const {
00400 scaleVec->push_back(1);
00401 for( int i=1; i<this->parNum_; ++i ) {
00402 scaleVec->push_back(0);
00403 }
00404 }
00405 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00406 double thisStep[] = {0.001, 0.01, 0.01, 0.1};
00407 TString thisParName[] = {"Pt scale", "Pt slope", "Phi ampl", "Phi phase"};
00408 double thisMini[] = {0.97, -0.1, -0.02, 0.};
00409 double thisMaxi[] = {1.03, 0.1, 0.02, 3.1416};
00410 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00411 }
00412 };
00413
00414
00415 template <class T>
00416 class scaleFunctionType12 : public scaleFunctionBase<T> {
00417 public:
00418 scaleFunctionType12() { this->parNum_ = 6; }
00419 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00420 return( (parScale[0] + parScale[1]*pt +
00421 parScale[2]*std::fabs(eta) +
00422 parScale[3]*eta*eta +
00423 (double)chg*parScale[4]*sin(phi+parScale[5]))*pt );
00424 }
00425
00426 virtual void resetParameters(std::vector<double> * scaleVec) const {
00427 scaleVec->push_back(1);
00428 for( int i=1; i<this->parNum_; ++i ) {
00429 scaleVec->push_back(0);
00430 }
00431 }
00432 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00433 double thisStep[] = {0.001, 0.01, 0.01, 0.01, 0.01, 0.1};
00434 TString thisParName[] = {"Pt scale", "Pt slope", "Eta slope", "Eta quadr", "Phi ampl", "Phi phase"};
00435 double thisMini[] = {0.97, -0.1, -0.2, -0.2, -0.02, 0.0};
00436 double thisMaxi[] = {1.03, 0.1, 0.2, 0.2, 0.02, 3.1416};
00437 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00438 }
00439 };
00440
00441
00442 template <class T>
00443 class scaleFunctionType13 : public scaleFunctionBase<T> {
00444 public:
00445 scaleFunctionType13() { this->parNum_ = 8; }
00446 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00447 if (chg>0) {
00448 return( (parScale[0] + parScale[1]*pt +
00449 parScale[2]*std::fabs(eta) +
00450 parScale[3]*eta*eta +
00451 parScale[4]*sin(phi+parScale[5]))*pt );
00452 }
00453
00454 return( (parScale[0] + parScale[1]*pt +
00455 parScale[2]*std::fabs(eta) +
00456 parScale[3]*eta*eta +
00457 parScale[6]*sin(phi+parScale[7]))*pt );
00458
00459 }
00460
00461 virtual void resetParameters(std::vector<double> * scaleVec) const {
00462 scaleVec->push_back(1);
00463 for( int i=1; i<this->parNum_; ++i ) {
00464 scaleVec->push_back(0);
00465 }
00466 }
00467 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00468 double thisStep[] = {0.001, 0.01, 0.01, 0.01, 0.01, 0.1, 0.01, 0.1};
00469 TString thisParName[] = {"Pt scale", "Pt slope", "Eta slope", "Eta quadr", "Phi ampl+", "Phi phase+", "Phi ampl-", "Phi phase-"};
00470 double thisMini[] = {0.99, -0.01, -0.02, -0.02, -0.02, 0.0, -0.02, 0.0};
00471 double thisMaxi[] = {1.01, 0.01, 0.02, 0.02, 0.02, 3.1416, 0.02, 3.1416};
00472 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00473 }
00474 };
00475
00476
00477
00478 template <class T>
00479 class scaleFunctionType14 : public scaleFunctionBase<T> {
00480 public:
00481 scaleFunctionType14() { this->parNum_ = 10; }
00482 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00483
00484
00485
00486
00487
00488
00489
00490
00491 return( ( parScale[0] +
00492 parScale[1]*pt + parScale[2]*pt*pt + parScale[3]*pt*pt*pt +
00493 parScale[4]*std::fabs(eta) + parScale[5]*eta*eta + parScale[6]*std::fabs(eta*eta*eta) +
00494 parScale[7]*eta*eta*eta*eta + parScale[8]*std::fabs(eta*eta*eta*eta*eta) +
00495 parScale[9]*eta*eta*eta*eta*eta*eta )*pt );
00496 }
00497
00498 virtual void resetParameters(std::vector<double> * scaleVec) const {
00499 scaleVec->push_back(1);
00500 for( int i=1; i<this->parNum_; ++i ) {
00501 scaleVec->push_back(0);
00502 }
00503 }
00504 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00505 double thisStep[] = {0.001,
00506 0.01, 0.01, 0.001,
00507 0.01, 0.00001, 0.0000001, 0.00000001, 0.00000001, 0.000000001};
00508 TString thisParName[] = {"Pt offset",
00509 "Pt slope", "Pt quadr", "Pt cubic",
00510 "Eta slope", "Eta quadr", "Eta cubic", "Eta quartic", "Eta fifth grade", "Eta sixth grade"};
00511 if( muonType == 1 ) {
00512 double thisMini[] = {0.9,
00513 -0.3, -0.3, -0.3,
00514 -0.3, -0.3, -0.01, -0.001, -0.0001, -0.00001};
00515 double thisMaxi[] = {1.1,
00516 0.3, 0.3, 0.3,
00517 0.3, 0.3, 0.01, 0.001, 0.0001, 0.00001};
00518 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00519 } else {
00520 double thisMini[] = {0.97,
00521 -0.1, -0.001, -0.001,
00522 -0.1, -0.1, -0.1, -0.0001, -0.00001, -0.000001};
00523 double thisMaxi[] = {1.03,
00524 0.1, 0.001, 0.001,
00525 0.1, 0.1, 0.1, 0.0001, 0.00001, 0.000001};
00526 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00527 }
00528 }
00529 };
00530
00531
00532
00533 template <class T>
00534 class scaleFunctionType15 : public scaleFunctionBase<T> {
00535 public:
00536 scaleFunctionType15() { this->parNum_ = 5; }
00537 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00538 if( pt > parScale[0] ) {
00539 return( ( parScale[1] + parScale[3]*std::fabs(eta) + parScale[4]*pow(eta,2) )*pt );
00540 }
00541 else {
00542 return( ( parScale[2] + parScale[3]*std::fabs(eta) + parScale[4]*pow(eta,2) )*pt );
00543 }
00544 }
00545
00546 virtual void resetParameters(std::vector<double> * scaleVec) const {
00547
00548 scaleVec->push_back(originalPtRegionSeparator_);
00549
00550 scaleVec->push_back(1);
00551 scaleVec->push_back(1);
00552 for( int i=1; i<this->parNum_; ++i ) {
00553 scaleVec->push_back(0);
00554 }
00555 }
00556 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00557 originalPtRegionSeparator_ = parScale[0];
00558 double thisStep[] = {0.001,
00559 0.01, 0.01,
00560 0.01, 0.00001};
00561 TString thisParName[] = {"Pt offset",
00562 "Pt slope 1", "Pt slope 2",
00563 "Eta slope", "Eta quadr"};
00564 if( muonType == 1 ) {
00565 double thisMini[] = {30.,
00566 0.9, 0.9,
00567 -0.3, -0.3};
00568 double thisMaxi[] = {50.,
00569 1.1, 1.1,
00570 0.3, 0.3};
00571 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00572 } else {
00573 double thisMini[] = {30.,
00574 0.97, 0.97,
00575 -0.1, -0.01};
00576 double thisMaxi[] = {50.,
00577 1.03, 1.03,
00578 0.1, 0.01};
00579 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00580 }
00581 }
00582 protected:
00583 double originalPtRegionSeparator_;
00584 };
00585
00586
00587
00588 template <class T>
00589 class scaleFunctionType16 : public scaleFunctionBase<T> {
00590 public:
00591 scaleFunctionType16() { this->parNum_ = 5; }
00592 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00593 return (parScale[0] + parScale[1]*std::fabs(eta)+ parScale[2]*eta*eta + parScale[3]*pt + parScale[4]/(pt*pt))*pt;
00594 }
00595
00596 virtual void resetParameters(std::vector<double> * scaleVec) const {
00597 scaleVec->push_back(1);
00598 for( int i=1; i<this->parNum_; ++i ) {
00599 scaleVec->push_back(0);
00600 }
00601 }
00602 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00603 double thisStep[] = {0.0000001, 0.00000001, 0.0000001, 0.00000001, 0.00000001};
00604 TString thisParName[] = {"offset", "linearEta", "parabEta", "linearPt", "coeffOverPt2"};
00605 if( muonType == 1 ) {
00606 double thisMini[] = {30.,
00607 0.9, 0.9,
00608 -0.3, -0.3};
00609 double thisMaxi[] = {50.,
00610 1.1, 1.1,
00611 0.3, 0.3};
00612 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00613 } else {
00614 double thisMini[] = {0.9, -0.000001, -0.000001, -0.00001, -0.00001};
00615 double thisMaxi[] = { 1.1, 0.01, 0.001, 0.001, 0.1};
00616 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00617 }
00618 }
00619 };
00620
00621 template <class T>
00622 class scaleFunctionType17 : public scaleFunctionBase<T> {
00623 public:
00624 scaleFunctionType17() { this->parNum_ = 4; }
00625 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00626 return (parScale[0]*std::fabs(eta)+ parScale[1]*eta*eta + pt/(parScale[2]*pt + parScale[3]))*pt;
00627 }
00628
00629
00630 virtual void resetParameters(std::vector<double> * scaleVec) const {
00631 scaleVec->push_back(1);
00632 for( int i=1; i<this->parNum_; ++i ) {
00633 scaleVec->push_back(0);
00634 }
00635 }
00636 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00637 double thisStep[] = {0.00000001, 0.000000001, 0.00000001, 0.00000001};
00638 TString thisParName[] = {"linearEta", "parabEta", "coeffPt", "coeffOverPt"};
00639 if( muonType == 1 ) {
00640 double thisMini[] = {30.,
00641 0.9, 0.9,
00642 -0.3};
00643 double thisMaxi[] = {50.,
00644 1.1, 1.1,
00645 0.3};
00646 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00647 } else {
00648 double thisMini[] = {-0.000001, -0.000001, 0.8, -0.001};
00649 double thisMaxi[] = { 0.01, 0.005, 1.2, 0.001};
00650 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00651 }
00652 }
00653 };
00654 template <class T>
00655 class scaleFunctionType18 : public scaleFunctionBase<T> {
00656 public:
00657 scaleFunctionType18() { this->parNum_ = 4; }
00658 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00659
00660 if(std::fabs(eta)<0.2)
00661 return parScale[0]*pt;
00662 else if(std::fabs(eta)>0.2 && std::fabs(eta)<1.1)
00663 return parScale[1]*pt;
00664 else if(std::fabs(eta)>1.1 && std::fabs(eta)<1.5)
00665 return parScale[2]*pt;
00666 else
00667 return parScale[3]*pt;
00668 }
00669
00670 virtual void resetParameters(std::vector<double> * scaleVec) const {
00671 for( int i=1; i<this->parNum_; ++i ) {
00672 scaleVec->push_back(1);
00673 }
00674 }
00675
00676 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00677 double thisStep[] = {0.00000001, 0.000000001, 0.00000001, 0.00000001};
00678 TString thisParName[] = {"etaCentr", "barrel", "overlap", "endcaps"};
00679 if( muonType == 1 ) {
00680 double thisMini[] = {30.,
00681 0.9, 0.9,
00682 -0.3};
00683 double thisMaxi[] = {50.,
00684 1.1, 1.1,
00685 0.3};
00686 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00687 } else {
00688 double thisMini[] = {0.9, 0.9, 0.9, 0.9};
00689 double thisMaxi[] = {1.1, 1.1, 1.1, 1.1};
00690 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00691 }
00692 }
00693 };
00694
00695
00696
00697
00698
00699
00700 template <class T>
00701 class scaleFunctionType19 : public scaleFunctionBase<T> {
00702 public:
00703 scaleFunctionType19() { this->parNum_ = 9; }
00704 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00705 if (chg>0) {
00706 return( (parScale[0] + parScale[1]*sin(parScale[2]*phi + parScale[3])+ parScale[4]*std::fabs(eta) + parScale[5]*eta*eta )*pt);
00707 }
00708 return( (parScale[0] + parScale[6]*sin(parScale[7]*phi + parScale[8])+ parScale[4]*std::fabs(eta) + parScale[5]*eta*eta )*pt );
00709 }
00710
00711
00712 virtual void resetParameters(std::vector<double> * scaleVec) const {
00713 scaleVec->push_back(1);
00714 for( int i=1; i<this->parNum_; ++i ) {
00715 scaleVec->push_back(0);
00716 }
00717 }
00718
00719 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const
00720 std::vector<int> & parScaleOrder, const int muonType)
00721 {
00722 double thisStep[] = {0.001, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
00723 TString thisParName[] = {"Phi offset", "Phi ampl Pos","Phi freq Pos", "Phi phase Pos","Eta slope", "Eta quadr","Phi ampl Neg","Phi freq Neg", "Phi phase Neg"};
00724 if( muonType == 1 ) {
00725 double thisMini[] = {0.9, -0.1, -0.1, -0.1, -0.02, -0.02, -0.1, -0.1, -0.1};
00726 double thisMaxi[] = {1.1, 0.1, 0.1, 0.1, 0.02, 0.02, 0.1, 0.1, 0.1};
00727 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00728 } else {
00729 double thisMini[] = {0.9, -0.1, -2.0, 0., -0.1, -0.01, -0.1, -2.0, 0. };
00730 double thisMaxi[] = {1.1, 0.1, 2.0, 6.28, 0.1, 0.01, 0.1, 2.0, 3.14 };
00731
00732 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00733 }
00734 }
00735 };
00736
00737
00738 template <class T>
00739 class scaleFunctionType20 : public scaleFunctionBase<T> {
00740 public:
00741 scaleFunctionType20() { this->parNum_ = 10; }
00742 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00743 if( pt < parScale[8] ) {
00744 return( (parScale[0] + parScale[3] + parScale[6]*std::fabs(eta) + parScale[7]*eta*eta )*pt);
00745 }
00746 else if( pt < parScale[9] ) {
00747 return( (parScale[1] + parScale[4] + parScale[6]*std::fabs(eta) + parScale[7]*eta*eta )*pt);
00748 }
00749 return( (parScale[2] + parScale[5] + parScale[6]*std::fabs(eta) + parScale[7]*eta*eta )*pt);
00750 }
00751
00752
00753 virtual void resetParameters(std::vector<double> * scaleVec) const {
00754 scaleVec->push_back(1);
00755 scaleVec->push_back(1);
00756 scaleVec->push_back(1);
00757 for( int i=1; i<this->parNum_-2; ++i ) {
00758 scaleVec->push_back(0);
00759 }
00760 scaleVec->push_back(this->originalTransition1_);
00761 scaleVec->push_back(this->originalTransition2_);
00762 }
00763
00764 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const
00765 std::vector<int> & parScaleOrder, const int muonType) {
00766
00767 originalTransition1_ = parScale[8];
00768 originalTransition2_ = parScale[9];
00769
00770 double thisStep[] = {0.001, 0.01, 0.01, 0.1, 0.01, 0.01, 0.001, 0.001, 0.1, 0.1};
00771
00772 TString thisParName[] = {"offset1", "offset2", "offset3", "linearPt1", "linearPt2", "linearPt3",
00773 "linearEta", "parabEta", "transition1", "transition2"};
00774 if( muonType == 1 ) {
00775 double thisMini[] = {0.9, 0.9, 0.9, -1., -1., -1., -1., -1., 0., 0.};
00776 double thisMaxi[] = {1.1, 1.1, 1.1, 1., 1., 1., 1., 1., 100., 100.};
00777 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00778 } else {
00779 double thisMini[] = {0.9, 0.9, 0.9, -1., -1., -1., -1., -1., 0., 0.};
00780 double thisMaxi[] = {1.1, 1.1, 1.1, 1., 1., 1., 1., 1., 100., 100.};
00781
00782 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00783 }
00784 }
00785
00786 protected:
00787
00788 double originalTransition1_;
00789 double originalTransition2_;
00790 };
00791
00792
00793
00794
00795 template <class T>
00796 class scaleFunctionType21 : public scaleFunctionBase<T> {
00797 public:
00798 scaleFunctionType21() { this->parNum_ = 8; }
00799 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
00800 double ptPart = parScale[0] + parScale[1]*pt;
00801 if( eta >= 0 ) {
00802 return( (ptPart+
00803 parScale[2]*eta +
00804 parScale[3]*eta*eta +
00805 parScale[4]*eta*eta*eta)*pt );
00806 }
00807 return( (ptPart +
00808 parScale[5]*(-eta) +
00809 parScale[6]*eta*eta +
00810 parScale[7]*(-eta*eta*eta))*pt );
00811 }
00812
00813 virtual void resetParameters(std::vector<double> * scaleVec) const {
00814 scaleVec->push_back(1);
00815 for( int i=1; i<this->parNum_; ++i ) {
00816 scaleVec->push_back(0);
00817 }
00818 }
00819 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
00820 double thisStep[] = {0.00001, 0.000001, 0.0000001, 0.0000001, 0.0000001, 0.0000001, 0.0000001, 0.0000001};
00821 TString thisParName[] = {"Pt offset", "Pt slope", "Eta slope pos eta", "Eta quadr pos eta", "Eta cubic pos eta", "Eta slope neg eta", "Eta quadr neg eta", "Eta cubic neg eta"};
00822 if( muonType == 1 ) {
00823 double thisMini[] = {0.9, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3};
00824 double thisMaxi[] = {1.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3};
00825 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00826 } else {
00827 double thisMini[] = {0.9, -0.002, -0.01, -0.005, -0.005, -0.01, -0.005, -0.005};
00828 double thisMaxi[] = {1.1, 0.002, 0.01, 0.005, 0.005, 0.01, 0.005, 0.005};
00829 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00830 }
00831 }
00832 };
00833
00834
00835
00836 template <class T>
00837 class scaleFunctionType22 : public scaleFunctionBase<T>
00838 {
00839 public:
00840 scaleFunctionType22() { this->parNum_ = 5; }
00841 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const
00842 {
00843
00844 int negChg = 0;
00845 if( chg > 0 ) {
00846 if( phi > 0 ) {
00847 return (parScale[0] + parScale[1]*TMath::Abs(phi)*sin(2*phi + parScale[2]))*pt;
00848 }
00849 else {
00850 return (parScale[0] + parScale[3]*TMath::Abs(phi)*sin(2*phi + parScale[4]))*pt;
00851 }
00852 }
00853 else if( chg < 0 ) {
00854 if( phi > 0 ) {
00855 return (parScale[0] - parScale[1+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[2+negChg]))*pt;
00856 }
00857 else {
00858 return (parScale[0] - parScale[3+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[4+negChg]))*pt;
00859 }
00860 }
00861 std::cout << "Error: we should not be here." << std::endl;
00862 exit(1);
00863 return 1;
00864 }
00865
00866 virtual void resetParameters(std::vector<double> * scaleVec) const
00867 {
00868 scaleVec->push_back(1);
00869 for( int i=1; i<this->parNum_; ++i ) {
00870 scaleVec->push_back(0);
00871 }
00872 }
00873 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType)
00874 {
00875 double thisStep[] = {0.0001, 0.0001, 0.01, 0.0001, 0.01};
00876
00877 TString thisParName[] = {"Phi offset",
00878 "amplitude pos phi", "phase pos phi",
00879 "amplitude neg phi", "phase neg phi"};
00880
00881
00882
00883
00884 if( muonType == 1 ) {
00885 double thisMini[] = {0.9, -0.3, -0.3, -0.3, -0.3};
00886
00887 double thisMaxi[] = {1.1, 0.3, 0.3, 0.3, 0.3};
00888
00889 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00890 } else {
00891 double thisMini[] = {0.9, -0.1, -3, -0.1, -3};
00892
00893 double thisMaxi[] = {1.1, 0.1, 3, 0.1, 3};
00894
00895 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00896 }
00897 }
00898 };
00899
00900
00901
00902
00903 template <class T>
00904 class scaleFunctionType23 : public scaleFunctionBase<T>
00905 {
00906 public:
00907 scaleFunctionType23() { this->parNum_ = 11; }
00908 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const
00909 {
00910
00911 int negChg = 4;
00912 if( chg > 0 ) {
00913 if( phi > 0 ) {
00914 return (parScale[0] + parScale[9]*etaCorrection(eta) + parScale[10]*eta*eta + parScale[1]*TMath::Abs(phi)*sin(2*phi + parScale[2]))*pt;
00915 }
00916 else {
00917 return (parScale[0] + parScale[9]*etaCorrection(eta) + parScale[10]*eta*eta + parScale[3]*TMath::Abs(phi)*sin(2*phi + parScale[4]))*pt;
00918 }
00919 }
00920 else if( chg < 0 ) {
00921 if( phi > 0 ) {
00922 return (parScale[0] + parScale[9]*etaCorrection(eta) + parScale[10]*eta*eta - parScale[1+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[2+negChg]))*pt;
00923 }
00924 else {
00925 return (parScale[0] + parScale[9]*etaCorrection(eta) + parScale[10]*eta*eta - parScale[3+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[4+negChg]))*pt;
00926 }
00927 }
00928 std::cout << "Error: we should not be here." << std::endl;
00929 exit(1);
00930 return 1;
00931 }
00932 double etaCorrection(const double & eta) const
00933 {
00934 double fabsEta = std::fabs(eta);
00935 if( fabsEta < 0.2) return -0.00063509;
00936 else if( fabsEta < 0.4 ) return -0.000585369;
00937 else if( fabsEta < 0.6 ) return -0.00077363;
00938 else if( fabsEta < 0.8 ) return -0.000547868;
00939 else if( fabsEta < 1.0 ) return -0.000954819;
00940 else if( fabsEta < 1.2 ) return -0.000162139;
00941 else if( fabsEta < 1.4 ) return 0.0026909;
00942 else if( fabsEta < 1.6 ) return 0.000684376;
00943 else if( fabsEta < 1.8 ) return -0.00174534;
00944 else if( fabsEta < 2.0 ) return -0.00177076;
00945 else if( fabsEta < 2.2 ) return 0.00117463;
00946 else if( fabsEta < 2.4 ) return 0.000985705;
00947 else if( fabsEta < 2.6 ) return 0.00163941;
00948 return 0.;
00949 }
00950
00951
00952 virtual void resetParameters(std::vector<double> * scaleVec) const
00953 {
00954 scaleVec->push_back(1);
00955 for( int i=1; i<this->parNum_; ++i ) {
00956 scaleVec->push_back(0);
00957 }
00958 }
00959 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType)
00960 {
00961 double thisStep[] = {0.0001,
00962 0.0001, 0.01, 0.0001, 0.01,
00963 0.0001, 0.01, 0.0001, 0.01,
00964 0.001,
00965 0.00001};
00966 TString thisParName[] = {"Phi offset",
00967
00968
00969 "amplitude pos charge pos phi", "phase pos charge pos phi",
00970 "amplitude pos charge neg phi", "phase pos charge neg phi",
00971 "amplitude neg charge pos phi", "phase neg charge pos phi",
00972 "amplitude neg charge neg phi", "phase neg charge neg phi",
00973 "amplitude of eta correction",
00974 "quadratic eta"};
00975 if( muonType == 1 ) {
00976 double thisMini[] = {0.9,
00977 -0.3, -0.3, -0.3, -0.3,
00978 -0.3, -0.3, -0.3, -0.3,
00979 -10.,
00980 -1.};
00981 double thisMaxi[] = {1.1,
00982 0.3, 0.3, 0.3, 0.3,
00983 0.3, 0.3, 0.3, 0.3,
00984 10.,
00985 1.};
00986 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00987 } else {
00988 double thisMini[] = {0.9,
00989 -0.1, -3, -0.1, -3,
00990 -0.1, -3, -0.1, -3,
00991 -10.,
00992 -1.};
00993 double thisMaxi[] = {1.1,
00994 0.1, 3, 0.1, 3,
00995 0.1, 3, 0.1, 3,
00996 10.,
00997 1.};
00998 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00999 }
01000 }
01001 };
01002
01003
01004
01005 template <class T>
01006 class scaleFunctionType24 : public scaleFunctionBase<T>
01007 {
01008 public:
01009 scaleFunctionType24() { this->parNum_ = 10; }
01010 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const
01011 {
01012
01013 int negChg = 4;
01014 if( chg > 0 ) {
01015 if( phi > 0 ) {
01016 return (parScale[0] + parScale[9]*etaCorrection(eta) + parScale[1]*TMath::Abs(phi)*sin(2*phi + parScale[2]))*pt;
01017 }
01018 else {
01019 return (parScale[0] + parScale[9]*etaCorrection(eta) + parScale[3]*TMath::Abs(phi)*sin(2*phi + parScale[4]))*pt;
01020 }
01021 }
01022 else if( chg < 0 ) {
01023 if( phi > 0 ) {
01024 return (parScale[0] + parScale[9]*etaCorrection(eta) - parScale[1+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[2+negChg]))*pt;
01025 }
01026 else {
01027 return (parScale[0] + parScale[9]*etaCorrection(eta) - parScale[3+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[4+negChg]))*pt;
01028 }
01029 }
01030 std::cout << "Error: we should not be here." << std::endl;
01031 exit(1);
01032 return 1;
01033 }
01034 double etaCorrection(const double & eta) const
01035 {
01036 if( eta < -2.6 ) return 0.;
01037 else if( eta < -2.4) return 0.00205594;
01038 else if( eta < -2.2) return 0.000880532;
01039 else if( eta < -2.0) return 0.0013714;
01040 else if( eta < -1.8) return -0.00153122;
01041 else if( eta < -1.6) return -0.000894437;
01042 else if( eta < -1.4) return 0.000883338;
01043 else if( eta < -1.2) return 0.0027599;
01044 else if( eta < -1.0) return 8.57009e-05;
01045 else if( eta < -0.8) return -0.00092294;
01046 else if( eta < -0.6) return -0.000492001;
01047 else if( eta < -0.4) return -0.000948406;
01048 else if( eta < -0.2) return -0.000478767;
01049 else if( eta < 0.0) return -0.0006909;
01050 else if( eta < 0.2) return -0.000579281;
01051 else if( eta < 0.4) return -0.000691971;
01052 else if( eta < 0.6) return -0.000598853;
01053 else if( eta < 0.8) return -0.000603736;
01054 else if( eta < 1.0) return -0.000986699;
01055 else if( eta < 1.2) return -0.00040998;
01056 else if( eta < 1.4) return 0.00262189;
01057 else if( eta < 1.6) return 0.000485414;
01058 else if( eta < 1.8) return -0.00259624;
01059 else if( eta < 2.0) return -0.00201031;
01060 else if( eta < 2.2) return 0.000977849;
01061 else if( eta < 2.5) return 0.00109088;
01062 else if( eta < 2.6) return 0.00122289;
01063 return 0.;
01064 }
01065
01066
01067 virtual void resetParameters(std::vector<double> * scaleVec) const
01068 {
01069 scaleVec->push_back(1);
01070 for( int i=1; i<this->parNum_; ++i ) {
01071 scaleVec->push_back(0);
01072 }
01073 }
01074 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType)
01075 {
01076 double thisStep[] = {0.0001,
01077 0.0001, 0.01, 0.0001, 0.01,
01078 0.0001, 0.01, 0.0001, 0.01,
01079 0.001};
01080 TString thisParName[] = {"Phi offset",
01081
01082
01083 "amplitude pos charge pos phi", "phase pos charge pos phi",
01084 "amplitude pos charge neg phi", "phase pos charge neg phi",
01085 "amplitude neg charge pos phi", "phase neg charge pos phi",
01086 "amplitude neg charge neg phi", "phase neg charge neg phi",
01087 "amplitude of eta correction"};
01088 if( muonType == 1 ) {
01089 double thisMini[] = {0.9,
01090 -0.3, -0.3, -0.3, -0.3,
01091 -0.3, -0.3, -0.3, -0.3,
01092 -10.};
01093 double thisMaxi[] = {1.1,
01094 0.3, 0.3, 0.3, 0.3,
01095 0.3, 0.3, 0.3, 0.3,
01096 10.};
01097 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01098 } else {
01099 double thisMini[] = {0.9,
01100 -0.1, -3, -0.1, -3,
01101 -0.1, -3, -0.1, -3,
01102 -10.};
01103 double thisMaxi[] = {1.1,
01104 0.1, 3, 0.1, 3,
01105 0.1, 3, 0.1, 3,
01106 10.};
01107 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01108 }
01109 }
01110 };
01111
01112
01113
01114
01115 template <class T>
01116 class scaleFunctionType25 : public scaleFunctionBase<T>
01117 {
01118 public:
01119 scaleFunctionType25() { this->parNum_ = 19; }
01120 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const
01121 {
01122
01123 int negChg = 4;
01124 if( chg > 0 ) {
01125 if( phi > 0 ) {
01126 return (parScale[0] + etaCorrection(eta, parScale) + parScale[1]*TMath::Abs(phi)*sin(2*phi + parScale[2]))*pt;
01127 }
01128 else {
01129 return (parScale[0] + etaCorrection(eta, parScale) + parScale[3]*TMath::Abs(phi)*sin(2*phi + parScale[4]))*pt;
01130 }
01131 }
01132 else if( chg < 0 ) {
01133 if( phi > 0 ) {
01134 return (parScale[0] + etaCorrection(eta, parScale) - parScale[1+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[2+negChg]))*pt;
01135 }
01136 else {
01137 return (parScale[0] + etaCorrection(eta, parScale) - parScale[3+negChg]*TMath::Abs(phi)*sin(2*phi + parScale[4+negChg]))*pt;
01138 }
01139 }
01140 std::cout << "Error: we should not be here." << std::endl;
01141 exit(1);
01142 return 1;
01143 }
01144 double etaCorrection(const double & eta, const T & parScale) const
01145 {
01146 if( eta < -2.06 ) return cubicEta(-2.06, parScale, 9);
01147 else if( eta < -1.06 ) return cubicEta(eta, parScale, 9);
01148 else if( eta < 1.1 ) return( parScale[13] + parScale[14]*eta );
01149 else if( eta < 2. ) return cubicEta(eta, parScale, 15);
01150 return cubicEta(2., parScale, 15);
01151 }
01152 double cubicEta(const double & eta, const T & parScale, const int shift ) const
01153 {
01154 return( parScale[shift] + parScale[shift+1]*eta + parScale[shift+2]*eta*eta + parScale[shift+3]*eta*eta*eta );
01155 }
01156
01157
01158 virtual void resetParameters(std::vector<double> * scaleVec) const
01159 {
01160 scaleVec->push_back(1);
01161 for( int i=1; i<this->parNum_; ++i ) {
01162 scaleVec->push_back(0);
01163 }
01164 }
01165 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType)
01166 {
01167 double thisStep[] = {0.0001,
01168 0.0001, 0.01, 0.0001, 0.01,
01169 0.0001, 0.01, 0.0001, 0.01,
01170 0.0001, 0.0001, 0.0001, 0.0001,
01171 0.0001, 0.0001,
01172 0.0001, 0.0001, 0.0001, 0.0001};
01173 TString thisParName[] = {"Phi offset",
01174
01175
01176 "amplitude pos charge pos phi", "phase pos charge pos phi",
01177 "amplitude pos charge neg phi", "phase pos charge neg phi",
01178 "amplitude neg charge pos phi", "phase neg charge pos phi",
01179 "amplitude neg charge neg phi", "phase neg charge neg phi",
01180 "etaNeg0", "etaNeg1", "etaNeg2", "etaNeg3",
01181 "etaCent0", "etaCent1",
01182 "etaPos0", "etaPos1", "etaPos2", "etaPos3"};
01183 if( muonType == 1 ) {
01184 double thisMini[] = {0.9,
01185 -0.3, -0.3, -0.3, -0.3,
01186 -0.3, -0.3, -0.3, -0.3,
01187 -0.1, -1., -1., -1.,
01188 -0.1, -1.,
01189 -0.1, -1., -1., -1.};
01190 double thisMaxi[] = {1.1,
01191 0.3, 0.3, 0.3, 0.3,
01192 0.3, 0.3, 0.3, 0.3,
01193 0.1, 1., 1., 1.,
01194 0.1, 1.,
01195 0.1, 1., 1., 1.};
01196 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01197 } else {
01198 double thisMini[] = {0.9,
01199 -0.1, -3, -0.1, -3,
01200 -0.1, -3, -0.1, -3,
01201 -0.1, -0.6, -0.5, -0.08,
01202 -0.1, -0.001,
01203 -0.1, -0.1, -0.4, -0.01};
01204 double thisMaxi[] = {1.1,
01205 0.1, 3, 0.1, 3,
01206 0.1, 3, 0.1, 3,
01207 0.1, 0.1, 0.1, 0.01,
01208 0.1, 0.002,
01209 0.1, 0.8, 0.1, 0.2};
01210 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01211 }
01212 }
01213 };
01214
01215
01216
01217 template <class T>
01218 class scaleFunctionType26 : public scaleFunctionBase<T> {
01219 public:
01220 scaleFunctionType26() { this->parNum_ = 9; }
01221 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
01222 double ptPart = parScale[0] + parScale[1]*pt;
01223 double fabsEta = std::fabs(eta);
01224
01225 if( fabsEta > parScale[8] ) {
01226 if( eta > 0 ) {
01227 return( (ptPart+
01228 parScale[2]*eta +
01229 parScale[3]*eta*eta)*pt );
01230 }
01231 else {
01232 return( (ptPart+
01233 parScale[4]*eta +
01234 parScale[5]*eta*eta)*pt );
01235 }
01236 }
01237 return( (ptPart + parScale[6]*fabsEta + parScale[7]*eta*eta)*pt );
01238 }
01239
01240 virtual void resetParameters(std::vector<double> * scaleVec) const {
01241 scaleVec->push_back(1);
01242 for( int i=1; i<this->parNum_; ++i ) {
01243 scaleVec->push_back(0);
01244 }
01245 }
01246 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
01247 double thisStep[] = {0.00001, 0.000001, 0.0000001, 0.0000001, 0.0000001, 0.0000001, 0.0000001, 0.00001, 0.000001};
01248 TString thisParName[] = {"Pt offset", "Pt slope", "Eta slope pos eta", "Eta quadr pos eta", "Eta slope neg eta", "Eta quadr neg eta", "Eta splope barrel", "Eta quadr barrel", "Eta corr region"};
01249 if( muonType == 1 ) {
01250 double thisMini[] = {0.9, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3, 0.3, 0.};
01251 double thisMaxi[] = {1.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3};
01252 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01253 } else {
01254 double thisMini[] = {0.9, -0.002, -0.01, -0.005, -0.005, -0.01, -0.01, -0.005, 0.};
01255 double thisMaxi[] = {1.1, 0.002, 0.01, 0.005, 0.005, 0.01, 0.01, 0.005, 2.4};
01256 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01257 }
01258 }
01259 };
01260
01261
01262
01263
01264
01265 template <class T>
01266 class scaleFunctionType27 : public scaleFunctionBase<T> {
01267 public:
01268 scaleFunctionType27() { this->parNum_ = 13; }
01269 virtual double scale(const double & pt, const double & eta, const double & phi, const int chg, const T & parScale) const {
01270 double ptPart = parScale[0] + parScale[1]*pt;
01271 double fabsEta = std::fabs(eta);
01272
01273 if( fabsEta > parScale[12] ) {
01274 if( eta > 0 ) {
01275 return( (ptPart+parScale[2]+
01276 parScale[3]*(fabsEta - parScale[5]) +
01277 parScale[4]*(fabsEta - parScale[5])*(fabsEta - parScale[5]))*pt );
01278 }
01279 else {
01280 return( (ptPart+parScale[6]+
01281 parScale[7]*(fabsEta - parScale[9]) +
01282 parScale[8]*(fabsEta - parScale[9])*(fabsEta - parScale[9]))*pt );
01283 }
01284 }
01285 return( (ptPart + parScale[10]*fabsEta + parScale[11]*eta*eta)*pt );
01286 }
01287
01288 virtual void resetParameters(std::vector<double> * scaleVec) const {
01289 scaleVec->push_back(1);
01290 for( int i=1; i<this->parNum_; ++i ) {
01291 scaleVec->push_back(0);
01292 }
01293 }
01294 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parScale, const std::vector<int> & parScaleOrder, const int muonType) {
01295 double thisStep[] = {0.00001, 0.000001,
01296 0.000001, 0.0000001, 0.0000001, 0.0000001,
01297 0.000001, 0.0000001, 0.0000001, 0.0000001,
01298 0.0000001, 0.0000001,
01299 0.00001};
01300 TString thisParName[] = {"Pt offset", "Pt slope",
01301 "Eta shift pos eta", "Eta slope pos eta", "Eta quadr pos eta", "Eta center pos eta",
01302 "Eta shift neg eta", "Eta slope neg eta", "Eta quadr neg eta", "Eta center neg eta",
01303 "Eta splope barrel", "Eta quadr barrel",
01304 "Eta corr region"};
01305 if( muonType == 1 ) {
01306 double thisMini[] = {0.9, -0.3,
01307 -0.3, -0.3, -0.3, -0.3,
01308 -0.3, -0.3, -0.3, -0.3,
01309 -0.3, -0.3,
01310 0.};
01311 double thisMaxi[] = {1.1, 0.3,
01312 0.3, 0.3, 0.3, 0.3,
01313 0.3, 0.3, 0.3, 0.3,
01314 0.3, 0.3,
01315 0.3};
01316 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01317 } else {
01318 double thisMini[] = {0.9, -0.002,
01319 -0.01, -0.01, -0.005, 0.,
01320 -0.01, -0.01, -0.005, 0.,
01321 -0.01, -0.005,
01322 0.};
01323 double thisMaxi[] = {1.1, 0.002,
01324 0.01, 0.01, 0.005, 2.4,
01325 0.01, 0.01, 0.005, 2.4,
01326 0.01, 0.005,
01327 2.4};
01328 this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
01329 }
01330 }
01331 };
01332
01333
01334
01336 scaleFunctionBase<double * > * scaleFunctionService( const int identifier );
01337
01339 scaleFunctionBase<std::vector<double> > * scaleFunctionVecService( const int identifier );
01340
01341
01342
01343
01344
01345 class smearFunctionBase {
01346 public:
01347 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) = 0;
01348 smearFunctionBase() {
01349 cotgth_ = 0.;
01350 gRandom_ = new TRandom();
01351 }
01352 virtual ~smearFunctionBase() = 0;
01353 protected:
01354 void smearEta(double & eta) {
01355 double theta;
01356 if (cotgth_!=0) {
01357 theta = atan(1/cotgth_);
01358 } else {
01359 theta = TMath::Pi()/2;
01360 }
01361 if (theta<0) theta += TMath::Pi();
01362 eta = -log(tan(theta/2));
01363 }
01364 double cotgth_;
01365 TRandom * gRandom_;
01366 };
01367 inline smearFunctionBase::~smearFunctionBase() { }
01368
01369
01370
01371 class smearFunctionType0 : public smearFunctionBase {
01372 public:
01373 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) { }
01374 };
01375
01376
01377 class smearFunctionType1 : public smearFunctionBase {
01378 public:
01379 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) {
01380 pt = pt*(1.0+y[0]*parSmear[0]*pt);
01381 phi = phi*(1.0+y[1]*parSmear[1]);
01382 double tmp = 2*atan(exp(-eta));
01383 cotgth_ = cos(tmp)/sin(tmp)*(1.0+y[2]*parSmear[2]);
01384 smearEta(eta);
01385 }
01386 };
01387
01388 class smearFunctionType2 : public smearFunctionBase {
01389 public:
01390 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) {
01391 pt = pt*(1.0+y[0]*parSmear[0]*pt+y[1]*parSmear[1]*std::fabs(eta));
01392 phi = phi*(1.0+y[2]*parSmear[2]);
01393 double tmp = 2*atan(exp(-eta));
01394 cotgth_ = cos(tmp)/sin(tmp)*(1.0+y[3]*parSmear[3]);
01395 smearEta(eta);
01396 }
01397 };
01398
01399 class smearFunctionType3 : public smearFunctionBase {
01400 public:
01401 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) {
01402 pt = pt*(1.0+y[0]*parSmear[0]*pt+y[1]*parSmear[1]*std::fabs(eta));
01403 phi = phi*(1.0+y[2]*parSmear[2]);
01404 double tmp = 2*atan(exp(-eta));
01405 cotgth_ = cos(tmp)/sin(tmp)*(1.0+y[3]*parSmear[3]+y[4]*parSmear[4]*std::fabs(eta));
01406 smearEta(eta);
01407 }
01408 };
01409
01410
01411
01412 class smearFunctionType4 : public smearFunctionBase {
01413 public:
01414 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) {
01415 pt = pt*(1.0+y[0]*parSmear[0]*pt+y[1]*parSmear[1]*std::fabs(eta)+y[5]*parSmear[5]*pow(pt,2));
01416 phi = phi*(1.0+y[2]*parSmear[2]);
01417 double tmp = 2*atan(exp(-eta));
01418 cotgth_ = cos(tmp)/sin(tmp)*(1.0+y[3]*parSmear[3]+y[4]*parSmear[4]*std::fabs(eta));
01419 smearEta(eta);
01420 }
01421 };
01422
01423 class smearFunctionType5 : public smearFunctionBase {
01424 public:
01425 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) {
01426 pt = pt*(1.0+y[0]*parSmear[0]*pt+y[1]*parSmear[1]*std::fabs(eta)+y[5]*parSmear[5]*pow(pt,2));
01427 phi = phi*(1.0+y[2]*parSmear[2]+y[6]*parSmear[6]*pt);
01428 double tmp = 2*atan(exp(-eta));
01429 cotgth_ = cos(tmp)/sin(tmp)*(1.0+y[3]*parSmear[3]+y[4]*parSmear[4]*std::fabs(eta));
01430 smearEta(eta);
01431 }
01432 };
01433
01434
01435 class smearFunctionType6 : public smearFunctionBase {
01436 public:
01437 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear) {
01438 double sigmaSmear = 0;
01439 double sigmaPtAl = 0;
01440 double sigmaPtMisal = 0;
01441 double ptPart = parSmear[0] + parSmear[1]*1/pt + pt*parSmear[2];
01442 double fabsEta = std::fabs(eta);
01443
01444 sigmaPtAl = parSmear[14]*etaByPoints(eta, parSmear[15]);
01445
01446 if (std::fabs(eta)<=1.4){
01447 sigmaPtMisal = ptPart + parSmear[3] + parSmear[4]*std::fabs(eta) + parSmear[5]*eta*eta;
01448 sigmaSmear = sqrt(std::fabs(pow(sigmaPtMisal,2)-pow(sigmaPtAl,2)));
01449 pt = pt*gRandom_->Gaus(1,sigmaSmear);
01450 }
01451 else if (eta>1.4){
01452 double par = parSmear[3] + parSmear[4]*1.4 + parSmear[5]*1.4*1.4 - (parSmear[6] + parSmear[7]*(1.4-parSmear[8]) + parSmear[9]*(1.4-parSmear[8])*(1.4-parSmear[8]));
01453 sigmaPtMisal = par + ptPart + parSmear[6] + parSmear[7]*std::fabs((fabsEta-parSmear[8])) + parSmear[9]*(fabsEta-parSmear[8])*(fabsEta-parSmear[8]);
01454 sigmaSmear = sqrt(std::fabs(pow(sigmaPtMisal,2)-pow(sigmaPtAl,2)));
01455 pt = pt*gRandom_->Gaus(1,sigmaSmear);
01456 }
01457 else{
01458 double par = parSmear[3] + parSmear[4]*1.4 + parSmear[5]*1.4*1.4 - (parSmear[10] + parSmear[11]*(1.4-parSmear[12]) + parSmear[13]*(1.4-parSmear[12])*(1.4-parSmear[12]));
01459 sigmaPtMisal = par + ptPart + parSmear[10] + parSmear[11]*std::fabs((fabsEta-parSmear[12])) + parSmear[13]*(fabsEta-parSmear[12])*(fabsEta-parSmear[12]);
01460 sigmaSmear = sqrt(std::fabs(pow(sigmaPtMisal,2)-pow(sigmaPtAl,2)));
01461 pt = pt*gRandom_->Gaus(1,sigmaSmear);
01462 }
01463 }
01464 protected:
01469 double etaByPoints(const double & inEta, const double & border) {
01470 Double_t eta = std::fabs(inEta);
01471 if( 0. <= eta && eta <= 0.2 ) return 0.00942984;
01472 else if( 0.2 < eta && eta <= 0.4 ) return 0.0104489;
01473 else if( 0.4 < eta && eta <= 0.6 ) return 0.0110521;
01474 else if( 0.6 < eta && eta <= 0.8 ) return 0.0117338;
01475 else if( 0.8 < eta && eta <= 1.0 ) return 0.0138142;
01476 else if( 1.0 < eta && eta <= 1.2 ) return 0.0165826;
01477 else if( 1.2 < eta && eta <= 1.4 ) return 0.0183663;
01478 else if( 1.4 < eta && eta <= 1.6 ) return 0.0169904;
01479 else if( 1.6 < eta && eta <= 1.8 ) return 0.0173289;
01480 else if( 1.8 < eta && eta <= 2.0 ) return 0.0205821;
01481 else if( 2.0 < eta && eta <= 2.2 ) return 0.0250032;
01482 else if( 2.2 < eta && eta <= 2.4 ) return 0.0339477;
01483
01484 else if( 2.4 < eta && eta <= 2.6 ) return border;
01485 return ( 0. );
01486 }
01487 };
01488
01489 class smearFunctionType7 : public smearFunctionBase
01490 {
01491 public:
01492 virtual void smear(double & pt, double & eta, double & phi, const double * y, const std::vector<double> & parSmear)
01493 {
01494 double sigmaSquared = sigmaPtDiff.squaredDiff(eta);
01495 TF1 G("G", "[0]*exp(-0.5*pow(x,2)/[1])", -5., 5.);
01496 double norm = 1/(sqrt(2*TMath::Pi()*sigmaSquared));
01497 G.SetParameter (0,norm);
01498 G.SetParameter (1,sigmaSquared);
01499 pt = pt*(1-G.GetRandom());
01500 }
01501 SigmaPtDiff sigmaPtDiff;
01502 };
01503
01505 smearFunctionBase * smearFunctionService( const int identifier );
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01521 template <class T>
01522 class resolutionFunctionBase {
01523 public:
01524 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) = 0;
01525 virtual double sigmaPtError(const double & pt, const double & eta, const T & parval, const T & parError)
01526 {
01527 return 0.;
01528 }
01529 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) = 0;
01530 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) = 0;
01531 virtual double covPt1Pt2(const double & pt1, const double & eta1, const double & pt2, const double & eta2, const T & parval)
01532 {
01533 return 0.;
01534 }
01535 resolutionFunctionBase() {}
01536 virtual ~resolutionFunctionBase() = 0;
01538 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) = 0;
01539 virtual int parNum() const { return parNum_; }
01540 protected:
01541 int parNum_;
01543 virtual void setPar(double* Start, double* Step, double* Mini, double* Maxi, int* ind,
01544 TString* parname, const T & parResol, const std::vector<int> & parResolOrder,
01545 double* thisStep, double* thisMini, double* thisMaxi, TString* thisParName ) {
01546 for( int iPar=0; iPar<this->parNum_; ++iPar ) {
01547 Start[iPar] = parResol[iPar];
01548 Step[iPar] = thisStep[iPar];
01549 Mini[iPar] = thisMini[iPar];
01550 Maxi[iPar] = thisMaxi[iPar];
01551 ind[iPar] = parResolOrder[iPar];
01552 parname[iPar] = thisParName[iPar];
01553 }
01554 }
01555 virtual void setPar(double* Start, double* Step, double* Mini, double* Maxi, int* ind,
01556 TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const std::vector<ParameterSet> & parSet ) {
01557 if( int(parSet.size()) != this->parNum_ ) {
01558 std::cout << "Error: wrong number of parameter initializations = " << parSet.size() << ". Number of parameters is " << this->parNum_ << std::endl;
01559 exit(1);
01560 }
01561 for( int iPar=0; iPar<this->parNum_; ++iPar ) {
01562 Start[iPar] = parResol[iPar];
01563 Step[iPar] = parSet[iPar].step;
01564 Mini[iPar] = parSet[iPar].mini;
01565 Maxi[iPar] = parSet[iPar].maxi;
01566 ind[iPar] = parResolOrder[iPar];
01567 parname[iPar] = parSet[iPar].name;
01568 }
01569 }
01570 };
01571 template <class T> inline resolutionFunctionBase<T>::~resolutionFunctionBase() { }
01572
01573
01574 template <class T>
01575 class resolutionFunctionType1 : public resolutionFunctionBase<T> {
01576 public:
01577 resolutionFunctionType1() { this->parNum_ = 3; }
01578 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
01579 return parval[0];
01580 }
01581 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
01582 return parval[1];
01583 }
01584 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01585 return parval[2];
01586 }
01587 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
01588 double thisStep[] = {0.001, 0.001, 0.001};
01589 TString thisParName[] = {"Pt res. sc.", "Phi res. sc.", "CotgThs res. sc."};
01590 double thisMini[] = {0., 0., 0.};
01591 if( muonType == 1 ) {
01592 double thisMaxi[] = {0.4, 0.4, 0.4};
01593 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01594 } else {
01595 double thisMaxi[] = {0.01, 0.02, 0.02};
01596 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01597 }
01598 }
01599 };
01600
01601
01602 template <class T>
01603 class resolutionFunctionType6 : public resolutionFunctionBase<T> {
01604 public:
01605 resolutionFunctionType6() { this->parNum_ = 15; }
01606 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
01607 return( parval[0]+parval[1]*pt+parval[2]*pow(pt,2)+parval[3]*pow(pt,3)+parval[4]*pow(pt,4)+parval[5]*std::fabs(eta)+parval[6]*pow(eta,2) );
01608 }
01609 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01610 return( parval[7]+parval[8]/pt+parval[9]*std::fabs(eta)+parval[10]*pow(eta,2) );
01611 }
01612 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
01613 return( parval[11]+parval[12]/pt+parval[13]*std::fabs(eta)+parval[14]*pow(eta,2) );
01614 }
01615 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
01616 double thisStep[] = { 0.005, 0.0005 ,0.000005 ,0.00000005 ,0.0000000005 ,0.0005 ,0.000005,
01617 0.000005 ,0.0005 ,0.00000005 ,0.0000005,
01618 0.00005 ,0.0005 ,0.00000005 ,0.000005};
01619 TString thisParName[] = { "Pt res. sc.", "Pt res. Pt sc.", "Pt res. Pt^2 sc.", "Pt res. Pt^3 sc.", "Pt res. Pt^4 sc.", "Pt res. Eta sc.", "Pt res. Eta^2 sc.",
01620 "Cth res. sc.", "Cth res. 1/Pt sc.", "Cth res. Eta sc.", "Cth res. Eta^2 sc.",
01621 "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc." };
01622 double thisMini[] = { 0.0, -0.01, -0.001, -0.01, -0.001, -0.001, -0.001,
01623 0.0, 0.0, -0.001, -0.001,
01624 0.0, 0.0, -0.001, -0.001 };
01625 if( muonType == 1 ) {
01626 double thisMaxi[] = { 1., 1., 1., 1., 1., 1., 1.,
01627 0.1, 1., 1., 1.,
01628 1., 1., 1., 1. };
01629 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01630 } else {
01631 double thisMaxi[] = { 0.1, 0.01, 0.01, 0.01, 0.01, 0.01, 0.1,
01632 0.01, 0.01, 0.01, 0.01,
01633 0.01, 0.01, 0.01, 0.01 };
01634 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01635 }
01636 }
01637 };
01638
01639
01640 template <class T>
01641 class resolutionFunctionType7 : public resolutionFunctionBase<T> {
01642 public:
01643 resolutionFunctionType7() { this->parNum_ = 12; }
01644
01645 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
01646 return( parval[0]+parval[1]*pt + parval[2]*std::fabs(eta)+parval[3]*pow(eta,2) );
01647 }
01648
01649 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01650 return( parval[4]+parval[5]/pt + parval[6]*std::fabs(eta)+parval[7]*pow(eta,2) );
01651 }
01652
01653 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
01654 return( parval[8]+parval[9]/pt + parval[10]*std::fabs(eta)+parval[11]*pow(eta,2) );
01655 }
01656 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
01657 double thisStep[] = { 0.002, 0.00002, 0.000002, 0.0002,
01658 0.00002, 0.0002, 0.0000002, 0.000002,
01659 0.00002, 0.0002, 0.00000002, 0.000002 };
01660 TString thisParName[] = { "Pt res. sc.", "Pt res. Pt sc.", "Pt res. Eta sc.", "Pt res. Eta^2 sc.",
01661 "Cth res. sc.", "Cth res. 1/Pt sc.", "Cth res. Eta sc.", "Cth res. Eta^2 sc.",
01662 "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc." };
01663 double thisMini[] = { 0.0, -0.01, -0.001, -0.0001,
01664 0.0, -0.001, -0.001, -0.00001,
01665 0.0, -0.001, -0.0001, -0.0001 };
01666 if( muonType == 1 ) {
01667 double thisMaxi[] = { 1., 1., 1., 1.,
01668 1., 1., 1., 0.1,
01669 1., 1., 1., 1. };
01670 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01671 } else {
01672 double thisMaxi[] = { 0.1, 0.01, 0.01, 0.1,
01673 0.01, 0.01, 0.1, 0.01,
01674 0.01, 0.01, 0.01, 0.01 };
01675 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01676 }
01677 }
01678 };
01679
01680
01681 template <class T>
01682 class resolutionFunctionType8 : public resolutionFunctionBase<T> {
01683 public:
01684 resolutionFunctionType8() { this->parNum_ = 12; }
01685
01686 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
01687 return( parval[0]+parval[1]*pt + parval[2]*etaByPoints(eta, parval[3]) );
01688 }
01689
01690 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01691 return( parval[4]+parval[5]/pt + parval[6]*std::fabs(eta)+parval[7]*eta*eta );
01692 }
01693
01694 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
01695 return( parval[8]+parval[9]/pt + parval[10]*std::fabs(eta)+parval[11]*eta*eta );
01696 }
01697 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
01698
01699 double thisStep[] = { 0.0000002, 0.0000001, 0.00001, 0.02,
01700 0.00002, 0.0002, 0.0000002, 0.00002,
01701 0.00002, 0.0002, 0.00000002, 0.000002 };
01702 TString thisParName[] = { "Pt res. sc.", "Pt res. Pt sc.", "Pt res. Eta sc.", "Pt res. eta border",
01703 "Cth res. sc.", "Cth res. 1/Pt sc.", "Cth res. Eta sc.", "Cth res. Eta^2 sc.",
01704 "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc." };
01705 double thisMini[] = { -0.03, -0.0000001, 0.1, 0.01,
01706 -0.001, 0.002, -0.0001, -0.0001,
01707 -0.0001, 0.0005, -0.0001, -0.00001 };
01708 if( muonType == 1 ) {
01709 double thisMaxi[] = { 1., 1., 1., 1.,
01710 1., 1., 1., 0.1,
01711 1., 1., 1., 1. };
01712 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01713 } else {
01714 double thisMaxi[] = { 0.03, 0.1, 1.4, 0.6,
01715 0.001, 0.005, 0.00004, 0.0007,
01716 0.001, 0.01, -0.0000015, 0.0004 };
01717 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01718 }
01719 }
01720 protected:
01725 double etaByPoints(const double & inEta, const double & border) {
01726 Double_t eta = std::fabs(inEta);
01727 if( 0. <= eta && eta <= 0.2 ) return 0.00942984;
01728 else if( 0.2 < eta && eta <= 0.4 ) return 0.0104489;
01729 else if( 0.4 < eta && eta <= 0.6 ) return 0.0110521;
01730 else if( 0.6 < eta && eta <= 0.8 ) return 0.0117338;
01731 else if( 0.8 < eta && eta <= 1.0 ) return 0.0138142;
01732 else if( 1.0 < eta && eta <= 1.2 ) return 0.0165826;
01733 else if( 1.2 < eta && eta <= 1.4 ) return 0.0183663;
01734 else if( 1.4 < eta && eta <= 1.6 ) return 0.0169904;
01735 else if( 1.6 < eta && eta <= 1.8 ) return 0.0173289;
01736 else if( 1.8 < eta && eta <= 2.0 ) return 0.0205821;
01737 else if( 2.0 < eta && eta <= 2.2 ) return 0.0250032;
01738 else if( 2.2 < eta && eta <= 2.4 ) return 0.0339477;
01739 else if( 2.4 < eta && eta <= 2.6 ) return border;
01740 return ( 0. );
01741 }
01742 };
01743
01744
01745 template <class T>
01746 class resolutionFunctionType9 : public resolutionFunctionBase<T> {
01747 public:
01748 resolutionFunctionType9() { this->parNum_ = 31; }
01749
01750 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
01751 double ptPart = 0.;
01752
01753 if( pt < 3 ) ptPart = parval[15];
01754 else if( pt < 4 ) ptPart = parval[16];
01755 else if( pt < 5 ) ptPart = parval[17];
01756 else if( pt < 6 ) ptPart = parval[18];
01757 else if( pt < 7 ) ptPart = parval[19];
01758 else if( pt < 8 ) ptPart = parval[20];
01759 else if( pt < 9 ) ptPart = parval[21];
01760 else if( pt < 10 ) ptPart = parval[22];
01761
01762 else ptPart = parval[0] + parval[1]*pt + parval[2]*pt*pt + parval[3]*pt*pt*pt + parval[4]*pt*pt*pt*pt;
01763
01764 double fabsEta = std::fabs(eta);
01765 double etaCoeff = parval[5];
01766 if( fabsEta < 0.1 ) etaCoeff = parval[23];
01767 else if( fabsEta < 0.2 ) etaCoeff = parval[24];
01768 else if( fabsEta < 0.3 ) etaCoeff = parval[25];
01769 else if( fabsEta < 0.4 ) etaCoeff = parval[26];
01770 else if( fabsEta < 0.5 ) etaCoeff = parval[27];
01771 else if( fabsEta < 0.6 ) etaCoeff = parval[28];
01772 else if( fabsEta < 0.7 ) etaCoeff = parval[29];
01773 else if( fabsEta < 0.8 ) etaCoeff = parval[30];
01774
01775 return( ptPart + etaCoeff*etaByPoints(eta, parval[6]) );
01776 }
01777
01778 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01779 return( parval[7]+parval[8]/pt + parval[9]*std::fabs(eta)+parval[10]*eta*eta );
01780 }
01781
01782 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
01783 return( parval[11]+parval[12]/pt + parval[13]*std::fabs(eta)+parval[14]*eta*eta );
01784 }
01785 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
01786
01787 double thisStep[] = { 0.0002, 0.000002, 0.0000002, 0.00000002, 0.000000002, 0.02, 0.02,
01788 0.00002, 0.0002, 0.0000002, 0.00002,
01789 0.00002, 0.0002, 0.00000002, 0.000002,
01790 0.001, 0.001, 0.001, 0.001,
01791 0.001, 0.001, 0.001, 0.001,
01792 0.001, 0.001, 0.001, 0.001,
01793 0.001, 0.001, 0.001, 0.001 };
01794 TString thisParName[] = { "Pt res. sc.", "Pt res. Pt sc.", "Pt res. Pt^2 sc.", "Pt res. Pt^3 sc.", "Pt res. Pt^4 sc",
01795 "Pt res. Eta sc.", "Pt res. eta border",
01796 "Cth res. sc.", "Cth res. 1/Pt sc.", "Cth res. Eta sc.", "Cth res. Eta^2 sc.",
01797 "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc.",
01798 "Pt by points sc. 0-3", "Pt by points sc. 3-4", "Pt by points sc. 4-5", "Pt by points sc. 5-6",
01799 "Pt by points sc. 6-7", "Pt by points sc. 7-8", "Pt by points sc. 8-9", "Pt by points sc. 9-10",
01800 "Eta scale for eta < 0.1", "Eta scale for eta < 0.2", "Eta scale for eta < 0.3", "Eta scale for eta < 0.4",
01801 "Eta scale for eta < 0.5", "Eta scale for eta < 0.6", "Eta scale for eta < 0.7", "Eta scale for eta < 0.8" };
01802 double thisMini[] = { -0.1, -0.001, -0.001, -0.001, -0.001, 0.001, 0.0001,
01803 -0.001, 0.002, -0.0001, -0.0001,
01804 -0.0001, 0.0005, -0.0001, -0.00001,
01805 -1., -1., -1., -1.,
01806 -1., -1., -1., -1.,
01807 -1., -1., -1., -1.,
01808 -1., -1., -1., -1. };
01809 if( muonType == 1 ) {
01810 double thisMaxi[] = { 1., 1., 1., 1., 1., 1., 1.,
01811 1., 1., 1., 0.1,
01812 1., 1., 1., 1.,
01813 1., 1., 1., 1.,
01814 1., 1., 1., 1.,
01815 1., 1., 1., 1.,
01816 3., 3., 3., 3. };
01817 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01818 } else {
01819 double thisMaxi[] = { 0.1, 0.001, 0.001, 0.001, 0.001, 1.5, 1.,
01820 0.001, 0.005, 0.00004, 0.0007,
01821 0.001, 0.01, -0.0000015, 0.0004,
01822 3., 3., 3., 3.,
01823 3., 3., 3., 3.,
01824 3., 3., 3., 3.,
01825 3., 3., 3., 3. };
01826 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01827 }
01828 }
01829 protected:
01834 double etaByPoints(const double & inEta, const double & border) {
01835 Double_t eta = std::fabs(inEta);
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879 if( 0. <= eta && eta <= 0.2 ) return 0.00942984;
01880 else if( 0.2 < eta && eta <= 0.4 ) return 0.0104489;
01881 else if( 0.4 < eta && eta <= 0.6 ) return 0.0110521;
01882 else if( 0.6 < eta && eta <= 0.8 ) return 0.0117338;
01883 else if( 0.8 < eta && eta <= 1.0 ) return 0.0138142;
01884 else if( 1.0 < eta && eta <= 1.2 ) return 0.0165826;
01885 else if( 1.2 < eta && eta <= 1.4 ) return 0.0183663;
01886 else if( 1.4 < eta && eta <= 1.6 ) return 0.0169904;
01887 else if( 1.6 < eta && eta <= 1.8 ) return 0.0173289;
01888 else if( 1.8 < eta && eta <= 2.0 ) return 0.0205821;
01889 else if( 2.0 < eta && eta <= 2.2 ) return 0.0250032;
01890 else if( 2.2 < eta && eta <= 2.4 ) return 0.0339477;
01891 else if( 2.4 < eta && eta <= 2.6 ) return border;
01892 return ( 0. );
01893
01894 }
01895 };
01896
01898
01899 template <class T>
01900 class resolutionFunctionType10 : public resolutionFunctionBase<T> {
01901 public:
01902 resolutionFunctionType10() { this->parNum_ = 21; }
01903
01904 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
01905 double fabsEta = std::fabs(eta);
01906 return( parval[0] + parval[1]*pt + parval[2]*pt*pt + parval[3]*pt*pt*pt + parval[4]*pt*pt*pt*pt
01907 + parval[5]*fabsEta + parval[6]*fabsEta*fabsEta + parval[7]*pow(fabsEta,3) + parval[8]*pow(fabsEta,4)
01908 + parval[9]*pow(fabsEta,5) + parval[10]*pow(fabsEta,6) + parval[11]*pow(fabsEta,7) + parval[12]*pow(fabsEta,8) );
01909 }
01910
01911 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01912 return( parval[13]+parval[14]/pt + parval[15]*std::fabs(eta)+parval[16]*eta*eta );
01913 }
01914
01915 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
01916 return( parval[17]+parval[18]/pt + parval[19]*std::fabs(eta)+parval[20]*eta*eta );
01917 }
01918 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
01919
01920 double thisStep[] = { 0.0002, 0.000002, 0.0000002, 0.00000002, 0.000000002,
01921 0.02, 0.02, 0.002, 0.0002,
01922 0.00002, 0.000002, 0.0000002, 0.00000002,
01923 0.00002, 0.0002, 0.0000002, 0.00002,
01924 0.00002, 0.0002, 0.00000002, 0.000002 };
01925 TString thisParName[] = { "Pt res. sc.", "Pt res. Pt sc.", "Pt res. Pt^2 sc.", "Pt res. Pt^3 sc.", "Pt res. Pt^4 sc",
01926 "Pt res. Eta sc.", "Pt res. Eta^2 sc." ,"Pt res. Eta^3 sc.", "Pt res. Eta^4 sc.",
01927 "Pt res. Eta^5 sc.", "Pt res. Eta^6 sc.", "Pt res. Eta^7 sc.", "Pt res. Eta^8 sc.",
01928 "Cth res. sc.", "Cth res. 1/Pt sc.", "Cth res. Eta sc.", "Cth res. Eta^2 sc.",
01929 "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc." };
01930 double thisMini[] = { -0.1, -0.001, -0.001, -0.001, -0.001,
01931 -2., -1., -0.1, -0.1,
01932 -0.1, -0.1, -0.1, -0.1,
01933 -0.001, 0.002, -0.0001, -0.0001,
01934 -0.0001, 0.0005, -0.0001, -0.00001,
01935 0.};
01936 if( muonType == 1 ) {
01937 double thisMaxi[] = { 1., 1., 1., 1., 1.,
01938 1., 1., 1., 1.,
01939 1., 1., 1., 1.,
01940 1., 1., 1., 0.1,
01941 1., 1., 1., 1. };
01942 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01943 } else {
01944 double thisMaxi[] = { 0.1, 0.001, 0.001, 0.001, 0.001,
01945 2., 1., 0.1, 0.1, 0.1,
01946 0.1, 0.1, 0.1, 0.1, 0.1,
01947 0.001, 0.005, 0.00004, 0.0007,
01948 0.001, 0.01, -0.0000015, 0.0004 };
01949 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01950 }
01951 }
01952 };
01953
01955
01956 template <class T>
01957 class resolutionFunctionType11 : public resolutionFunctionBase<T> {
01958 public:
01959 resolutionFunctionType11() { this->parNum_ = 8; }
01960
01961 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
01962 double fabsEta = std::fabs(eta);
01963 if(fabsEta<1.2)
01964 return (parval[0]+ parval[2]*1./pt + pt/(pt+parval[3]) + parval[4]*fabsEta + parval[5]*eta*eta);
01965 else
01966 return (parval[1]+ parval[2]*1./pt + pt/(pt+parval[3]) + parval[6]*std::fabs((fabsEta-1.6)) + parval[7]*(fabsEta-1.6)*(fabsEta-1.6));
01967 }
01968
01969 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01970 return( 0.004 );
01971 }
01972
01973 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
01974 return( 0.001 );
01975 }
01976 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
01977
01978 double thisStep[] = { 0.00001, 0.00001, 0.0000001, 0.00000001,
01979 0.00000001, 0.00000001, 0.00000001, 0.00000001 };
01980 TString thisParName[] = { "offsetEtaCentral", "offsetEtaHigh", "coeffOverPt", "coeffHighPt", "linaerEtaCentral", "parabEtaCentral", "linaerEtaHigh", "parabEtaHigh" };
01981 double thisMini[] = { -1.1, -1.1, -0.1, -0.1 , 0.0001, 0.0005, 0.0005, 0.001};
01982 if( muonType == 1 ) {
01983 double thisMaxi[] = { 1., 1., 1., 1.,
01984 1., 1., 1., 1.};
01985 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01986 } else {
01987 double thisMaxi[] = { -0.8, -0.8, -0.001, -0.001 , 0.005, 0.05, 0.05, 0.05};
01988 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
01989 }
01990 }
01991 };
01992
02000
02001 template <class T>
02002 class resolutionFunctionType12 : public resolutionFunctionBase<T> {
02003 public:
02004 resolutionFunctionType12() { this->parNum_ = 15; }
02005
02006 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02007 double fabsEta = std::fabs(eta);
02008
02009 double ptPart = parval[2]*1./pt + pt/(pt+parval[3]) + pt*parval[9] + pt*pt*parval[10];
02010
02011 if(fabsEta<parval[0]) {
02012
02013 double par = parval[1] + parval[6]*std::fabs((parval[0]-parval[8])) + parval[7]*(parval[0]-parval[8])*(parval[0]-parval[8]) - (parval[4]*parval[0] + parval[5]*parval[0]*parval[0]);
02014 return( par + ptPart + parval[4]*fabsEta + parval[5]*eta*eta );
02015 }
02016 else {
02017 return( parval[1]+ ptPart + parval[6]*std::fabs((fabsEta-parval[8])) + parval[7]*(fabsEta-parval[8])*(fabsEta-parval[8]) );
02018 }
02019 }
02020
02021
02022 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02023 return( parval[11]+parval[12]/pt + parval[13]*std::fabs(eta)+parval[14]*eta*eta );
02024 }
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02038 return( 0.001 );
02039 }
02040
02041 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
02042
02043 double thisStep[] = { 0.001, 0.00001, 0.0000001, 0.00000001,
02044 0.00000001, 0.00000001, 0.00000001, 0.00000001,
02045 0.001, 0.0001, 0.000001,
02046 0.00002, 0.0002, 0.0000002, 0.00002 };
02047
02048 TString thisParName[] = { "etaTransition", "offsetEtaHigh", "coeffOverPt", "coeffHighPt",
02049 "linaerEtaCentral", "parabEtaCentral", "linaerEtaHigh", "parabEtaHigh",
02050 "secondParabolaCenter", "linearPt", "quadraticPt",
02051 "Cth res. sc.", "Cth res. 1/Pt sc.", "Cth res. Eta sc.", "Cth res. Eta^2 sc." };
02052
02053 double thisMini[] = { 0.8, -1.1, 0., -1.1,
02054 0., 0.0005, 0.0005, 0.001,
02055 1.4, 0., 0.,
02056
02057 -0.1, 0., -0.1, -0.1 };
02058
02059 if( muonType == 1 ) {
02060 double thisMaxi[] = { 1., 1., 1., 1.,
02061 1., 1., 1., 1.,
02062 1., 1., 1.,
02063 1., 1., 1., 1. };
02064
02065
02066 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02067 } else {
02068 double thisMaxi[] = { 1.8, 0.8, 0.1, 0.1,
02069 0.005, 0.05, 0.05, 0.05,
02070 2.4, 2.0, 2.0,
02071
02072 0.1, 0.05, 0.1, 0.1 };
02073
02074 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02075 }
02076 }
02077 };
02078
02084
02085 template <class T>
02086 class resolutionFunctionType13 : public resolutionFunctionBase<T> {
02087 public:
02088 resolutionFunctionType13() { this->parNum_ = 15; }
02089
02090 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02091 double fabsEta = std::fabs(eta);
02092
02093 double ptPart = parval[2]*1./pt + pt/(pt+parval[3]) + pt*parval[13] + pt*pt*parval[14];
02094
02095 if(fabsEta<parval[9]) {
02096
02097 double par2 = parval[1] + parval[6]*std::fabs(parval[9] - parval[8]) + parval[7]*pow(parval[9] - parval[8], 2) - parval[10]*std::fabs(parval[9] - parval[12]) - parval[11]*pow(parval[9] - parval[12], 2);
02098 if( fabsEta<parval[0]) {
02099 double par1 = par2 + parval[10]*std::fabs(parval[0] - parval[12]) + parval[11]*pow(parval[0] - parval[12], 2) - parval[4]*parval[0] - parval[5]*parval[0]*parval[0];
02100 return( par1 + ptPart + parval[4]*fabsEta + parval[5]*eta*eta );
02101
02102 }
02103
02104 return( par2+ ptPart + parval[10]*std::fabs((fabsEta-parval[12])) + parval[11]*(fabsEta-parval[12])*(fabsEta-parval[12]) );
02105 }
02106 return( parval[1]+ ptPart + parval[6]*std::fabs((fabsEta-parval[8])) + parval[7]*(fabsEta-parval[8])*(fabsEta-parval[8]) );
02107 }
02108
02109 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02110 return( 0.004 );
02111 }
02112
02113 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02114 return( 0.001 );
02115 }
02116 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
02117
02118 double thisStep[] = { 0.001, 0.00001, 0.0000001, 0.00000001,
02119 0.00000001, 0.00000001, 0.00000001, 0.00000001,
02120 0.00000001, 0.00000001, 0.00000001, 0.00000001,
02121 0.001, 0.0001, 0.000001 };
02122
02123 TString thisParName[] = { "etaTransition1", "offsetEtaHigh", "coeffOverPt", "coeffHighPt",
02124 "linearEtaCentral", "parabEtaCentral", "linearEtaHighEta", "parabEtaHighEta",
02125 "centerEtaHighEta", "etaTransition2", "linearEtaSecondEta", "parabEtaSecondEta",
02126 "centerEtaSecondEta", "linearPt", "quadraticPt" };
02127
02128 double thisMini[] = { 0.8, -1.1, -1.1, -1.1,
02129 -1., 0.0005, 0.0005, 0.001,
02130 0.8, 1.1, 0.0005, 0.001,
02131 0.0, -1.0, -1.0 };
02132
02133 if( muonType == 1 ) {
02134 double thisMaxi[] = { 1., 1., 1., 1.,
02135 1., 1., 1., 1.,
02136 1., 1., 1., 1.,
02137 1., 1., 1. };
02138
02139
02140 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02141 } else {
02142 double thisMaxi[] = { 1.0, -0.8, 0.1, 0.1,
02143 0.005, 0.05, 0.05, 0.05,
02144 2.4, 1.4, 0.05, 0.05,
02145 1.8, 2.0, 2.0 };
02146
02147 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02148 }
02149 }
02150 };
02151
02152 template <class T>
02153 class resolutionFunctionType14 : public resolutionFunctionBase<T> {
02154 public:
02155 resolutionFunctionType14() { this->parNum_ = 6; }
02156
02157 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02158 double fabsEta = std::fabs(eta);
02159 if(fabsEta<1.2)
02160 return (parval[0] + parval[2]*fabsEta + parval[3]*eta*eta);
02161 else
02162 return (parval[1]+ parval[4]*std::fabs((fabsEta-1.6)) + parval[5]*(fabsEta-1.6)*(fabsEta-1.6));
02163 }
02164
02165 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02166 return( 0.004 );
02167 }
02168
02169 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02170 return( 0.001 );
02171 }
02172 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
02173
02174 double thisStep[] = { 0.00001, 0.00001, 0.0000001, 0.00000001, 0.00000001, 0.00000001, 0.00000001 };
02175 TString thisParName[] = { "offsetEtaCentral", "offsetEtaHigh", "linaerEtaCentral", "parabEtaCentral", "linaerEtaHigh", "parabEtaHigh" };
02176 double thisMini[] = { 0.005, 0.012, -0.01, 0., -0.05, 0.};
02177 if( muonType == 1 ) {
02178 double thisMaxi[] = { 1., 1., 1., 1., 1.,
02179 1.};
02180 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02181 } else {
02182 double thisMaxi[] = { 0.01, 0.02, 0.01, 0.05, 0.05, 0.1};
02183 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02184 }
02185 }
02186 protected:
02187 };
02188
02189
02190
02191 template <class T>
02192 class resolutionFunctionType15 : public resolutionFunctionBase<T> {
02193 public:
02194 resolutionFunctionType15() { this->parNum_ = 7; }
02195
02196 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02197 double fabsEta = std::fabs(eta);
02198 double ptPart = pt*parval[0];
02199
02200 if(fabsEta<=0.6)
02201 return (ptPart);
02202 else if(fabsEta>0.6 && fabsEta<=1.3) {
02203 double par = - parval[1]*0.6*0.6;
02204 return( par + ptPart + parval[1]*eta*eta );
02205 }
02206 else if (eta>1.3){
02207 double par = parval[1]*1.3*1.3 - (parval[2]*(1.3-parval[3])*(1.3-parval[3]));
02208 return( par + ptPart + parval[2]*(fabsEta-parval[3])*(fabsEta-parval[3]) );
02209 }
02210 else{
02211 double par = parval[1]*1.3*1.3 - (parval[4]*(1.3-parval[5]) + parval[6]*(1.3-parval[5])*(1.3-parval[5]));
02212 return( par + ptPart + parval[4]*std::fabs((fabsEta-parval[5])) + parval[6]*(fabsEta-parval[5])*(fabsEta-parval[5]) );
02213 }
02214 }
02215
02216 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02217 if (std::fabs(eta)<=1.2) return(0.000949148 );
02218 else if(eta<-1.2) return(-0.00645458 + -0.00579458*eta);
02219 else return(-0.00306283 + 0.00346136*eta);
02220 }
02221
02222 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02223 return( 0.000211 + 0.00001*std::fabs(eta) + 0.0000789*eta*eta);
02224 }
02225 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
02226 double thisStep[] = { 0.000001, 0.0000001,
02227 0.000001, 0.001,
02228 0.00000001, 0.001, 0.0001};
02229 TString thisParName[] = { "linearPt", "parabEtaCentral",
02230 "parabolicEtaRight", "rightParabCenter",
02231 "linearEtaLeft", "leftParabCenter", "parabolicEtaLeft" };
02232 double thisMini[] = { 0.00001, 0.0001,
02233 0.005, -5,
02234 -0.00006, 0.1, 0.002
02235 };
02236
02237 if( muonType == 1 ) {
02238 double thisMaxi[] = { 1., 1.,
02239 1., 1.,
02240 1., 1. ,1.,
02241 };
02242 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02243 } else {
02244 double thisMaxi[] = { 0.0005, 0.05,
02245 0.15, 1.99,
02246 0.005, 1.99, 0.15,
02247 };
02248
02249 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02250 }
02251 }
02252 };
02253
02254
02255 template <class T>
02256 class resolutionFunctionType17 : public resolutionFunctionBase<T> {
02257 public:
02258 resolutionFunctionType17() { this->parNum_ = 18; }
02259
02260 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02261 double fabsEta = std::fabs(eta);
02262 double ptPartBar = parval[0] + pt*parval[2];
02263 double ptPartOvlap = parval[16] + pt*parval[17];
02264 double ptPartEndc = parval[14] + pt*parval[15];
02265 if(fabsEta<=0.9) {
02266 return( ptPartBar + parval[3] + parval[4]*fabsEta);
02267 }
02268 else if( (eta > 0.9 && eta <= 1.4) || (eta < -0.9 && eta > -1.4)){
02269 return( ptPartOvlap + parval[3] + parval[4]*eta + parval[5]*eta*eta);
02270 }
02271 else if (eta>1.4){
02272 double par = parval[3] + parval[4]*1.4 + parval[5]*1.4*1.4 - (parval[6] + parval[7]*(1.4-parval[8]) + parval[9]*(1.4-parval[8])*(1.4-parval[8]));
02273 return( par + ptPartEndc + parval[6] + parval[7]*std::fabs((fabsEta-parval[8])) + parval[9]*(fabsEta-parval[8])*(fabsEta-parval[8]) );
02274 }
02275 else {
02276 double par = parval[3] + parval[4]*1.4 + parval[5]*1.4*1.4 - (parval[10] + parval[11]*(1.4-parval[12]) + parval[13]*(1.4-parval[12])*(1.4-parval[12]));
02277 return( par + ptPartEndc + parval[10] + parval[11]*std::fabs((fabsEta-parval[12])) + parval[13]*(fabsEta-parval[12])*(fabsEta-parval[12]) );
02278 }
02279 }
02280
02281 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02282 return( 0.004 );
02283 }
02284
02285 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02286 return( 0.001 );
02287 }
02288 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
02289 double thisStep[] = { 0.0001, 0.00001, 0.00001,
02290 0.001, 0.0001, 0.0000001,
02291 0.01, 0.001, 0.01, 0.001,
02292 0.01, 0.001, 0.01, 0.001,
02293 0.01, 0.00001, 0.01, 0.00001};
02294 TString thisParName[] = { "offsetPt", "hyperbolicPt", "linearPt",
02295 "offsetEtaCentral", "linaerEtaCentral", "parabEtaCentral",
02296 "offsetEtaEndcapRight", "linearEtaRight", "rightParabCenter", "parabolicEtaRight",
02297 "offsetEtaEndcapLeft", "linearEtaLeft", "leftParabCenter", "parabolicEtaLeft",
02298 "offsetPtEndc", "linearPtEndc", "offsetPtOvlap", "linearPtOvlap"
02299 };
02300 double thisMini[] = { -0.15, -0.001, 0.00005,
02301 -0.05, -0.1, 0.0,
02302 -0.6, -0.0009, 0., 0.0005,
02303 -0.6, -0.1, 1., 0.01,
02304 -1.5, 0.00005, -1.5, 0.00005
02305 };
02306
02307 if( muonType == 1 ) {
02308 double thisMaxi[] = { 1., 1., 1.,
02309 1., 1., 1.,
02310 1., 1. ,1., 1.,
02311 1., 1., 1., 1.,
02312 1., 1., 1., 1.,
02313 };
02314 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02315 } else {
02316 double thisMaxi[] = { 0.15, 0.8, 0.005,
02317 0.05, 0.1, 0.08,
02318 0.9, 0.5, 1.99, 0.15,
02319 0.9, 0.5, 1.99, 0.15,
02320 1.1, 0.005, 1.1, 0.005
02321 };
02322
02323 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02324 }
02325 }
02326 };
02327
02328 template <class T>
02329 class resolutionFunctionType18 : public resolutionFunctionBase<T> {
02330 public:
02331 resolutionFunctionType18() { this->parNum_ = 14; }
02332
02333 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02334 double fabsEta = std::fabs(eta);
02335 double ptPart = parval[0] + parval[1]*1./pt + pt*parval[2];
02336
02337 if(fabsEta<=0.6)
02338 return( ptPart + parval[3]);
02339 else if((eta>0.6 && eta<=1.3) || (eta>=-1.3 && eta<-0.6)) {
02340 double par = parval[3] - 0.6*parval[4] - 0.6*0.6*parval[5];
02341 return( ptPart + par + parval[4]*fabsEta + parval[5]*eta*eta );
02342 }
02343 else if (eta>1.3){
02344 double par = parval[3] - 0.6*parval[4] - 0.6*0.6*parval[5] + parval[4]*1.3 + parval[5]*1.3*1.3 - (parval[6] + parval[7]*std::fabs(1.3-parval[8]) + parval[9]*(1.3-parval[8])*(1.3-parval[8]));
02345 return( par + ptPart + parval[6] + parval[7]*std::fabs((fabsEta-parval[8])) + parval[9]*(fabsEta-parval[8])*(fabsEta-parval[8]) );
02346 }
02347 else{
02348 double par = parval[3] - 0.6*parval[4] - 0.6*0.6*parval[5] + parval[4]*1.3 + parval[5]*1.3*1.3 - (parval[10] + parval[11]*std::fabs(1.3-parval[12]) + parval[13]*(1.3-parval[12])*(1.3-parval[12]));
02349 return( par + ptPart + parval[10] + parval[11]*std::fabs((fabsEta-parval[12])) + parval[13]*(fabsEta-parval[12])*(fabsEta-parval[12]) );
02350 }
02351 }
02352
02353 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02354 return( 0.004 );
02355 }
02356
02357 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02358 return( 0.001 );
02359 }
02360 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
02361 double thisStep[] = { 0.01, 0.0001, 0.00001,
02362 0.001, 0.0001, 0.000001,
02363 0.01, 0.001, 0.01, 0.001,
02364 0.01, 0.001, 0.01, 0.001};
02365 TString thisParName[] = { "offsetPt", "hyperbolicPt", "linearPt",
02366 "offsetEtaCentral", "linaerEtaCentral", "parabEtaCentral",
02367 "offsetEtaEndcapRight", "linearEtaRight", "rightParabCenter", "parabolicEtaRight",
02368 "offsetEtaEndcapLeft", "linearEtaLeft", "leftParabCenter", "parabolicEtaLeft" };
02369 double thisMini[] = { -1.5, -0.001, 0.00005,
02370 -0.05, -0.1, 0.0,
02371 -0.6, -0.0009, 0., 0.0005,
02372 -0.6, -0.1, 1., 0.01
02373 };
02374
02375 if( muonType == 1 ) {
02376 double thisMaxi[] = { 1., 1., 1.,
02377 1., 1., 1.,
02378 1., 1. ,1., 1.,
02379 1., 1., 1., 1.
02380 };
02381 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02382 } else {
02383 double thisMaxi[] = { 1.1, 0.8, 0.005,
02384 0.05, 0.1, 0.08,
02385 0.9, 0.5, 1.99, 0.15,
02386 0.9, 0.5, 1.99, 0.15
02387 };
02388
02389 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02390 }
02391 }
02392 };
02393
02394
02395
02396
02397 template <class T>
02398 class resolutionFunctionType19 : public resolutionFunctionBase<T> {
02399 public:
02400 resolutionFunctionType19() { this->parNum_ = 4; }
02401
02402 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02403 double value = parval[0]+parval[1]*pt + parval[2]*etaByPoints(eta, parval[3]);
02404 if( value != value ) {
02405 std::cout << "parval[0] = " << parval[0] << ", parval[1]*"<<pt<<" = " << parval[1]*pt << "parval[2] = " << parval[1] << ",etaByPoints("<<eta<<", "<<parval[3]<<") = " << etaByPoints(eta, parval[3]) << std::endl;
02406 }
02407 return( parval[0] + parval[1]*pt + parval[2]*etaByPoints(eta, parval[3]) );
02408 }
02409
02410 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02411 return( 0.00043 + 0.0041/pt + (2.8e-06)*std::fabs(eta) + (7.7e-05)*eta*eta );
02412 }
02413
02414 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02415 return( 0.00011 + 0.0018/pt - (9.4e-07)*std::fabs(eta) + (2.2e-05)*eta*eta );
02416 }
02417 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType)
02418 {
02419 double thisStep[] = { 0.0000002, 0.0000001, 0.00001, 0.001 };
02420 TString thisParName[] = { "Pt res. sc.", "Pt res. Pt sc.", "Pt res. Eta sc.", "Pt res. eta border"};
02421 double thisMini[] = { -0.03, -0.0000001, 0.001, 0.01};
02422 if( muonType == 1 ) {
02423 double thisMaxi[] = { 1., 1., 1., 1.};
02424 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02425 } else {
02426 double thisMaxi[] = { 0.03, 0.1, 2., 0.6};
02427 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02428 }
02429 }
02430 protected:
02435 double etaByPoints(const double & inEta, const double & border) {
02436 Double_t eta = std::fabs(inEta);
02437 if( 0. <= eta && eta <= 0.2 ) return 0.00942984;
02438 else if( 0.2 < eta && eta <= 0.4 ) return 0.0104489;
02439 else if( 0.4 < eta && eta <= 0.6 ) return 0.0110521;
02440 else if( 0.6 < eta && eta <= 0.8 ) return 0.0117338;
02441 else if( 0.8 < eta && eta <= 1.0 ) return 0.0138142;
02442 else if( 1.0 < eta && eta <= 1.2 ) return 0.0165826;
02443 else if( 1.2 < eta && eta <= 1.4 ) return 0.0183663;
02444 else if( 1.4 < eta && eta <= 1.6 ) return 0.0169904;
02445 else if( 1.6 < eta && eta <= 1.8 ) return 0.0173289;
02446 else if( 1.8 < eta && eta <= 2.0 ) return 0.0205821;
02447 else if( 2.0 < eta && eta <= 2.2 ) return 0.0250032;
02448 else if( 2.2 < eta && eta <= 2.4 ) return 0.0339477;
02449 else if( 2.4 < eta && eta <= 2.6 ) return border;
02450 return ( 0. );
02451 }
02452 };
02453
02454 template <class T>
02455 class resolutionFunctionType20 : public resolutionFunctionBase<T> {
02456 public:
02457 resolutionFunctionType20() { this->parNum_ = 9; }
02458
02459 virtual double sigmaPt(const double & pt, const double & eta, const T & parval) {
02460 double fabsEta = std::fabs(eta);
02461
02462 if(fabsEta<parval[0]) {
02463
02464 double par = parval[1] + parval[4]*std::fabs((parval[0]-parval[6])) + parval[5]*(parval[0]-parval[6])*(parval[0]-parval[6]) - (parval[2]*parval[0] + parval[3]*parval[0]*parval[0]);
02465 return( par + parval[2]*fabsEta + parval[3]*eta*eta );
02466 }
02467 else {
02468 return( parval[1]+ parval[4]*std::fabs((fabsEta-parval[6])) + parval[5]*(fabsEta-parval[6])*(fabsEta-parval[6]) );
02469 }
02470 }
02471
02472
02473 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02474 return( parval[7]+parval[8]/pt );
02475 }
02476
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02489 return( 0.001 );
02490 }
02491
02492 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType) {
02493
02494 double thisStep[] = { 0.001, 0.00001,
02495 0.00000001, 0.00000001, 0.00000001, 0.00000001,
02496 0.001,
02497 0.00002, 0.0002 };
02498
02499 TString thisParName[] = { "etaTransition", "offsetEtaHigh",
02500 "linaerEtaCentral", "parabEtaCentral", "linaerEtaHigh", "parabEtaHigh",
02501 "secondParabolaCenter",
02502 "Cth res. sc.", "Cth res. 1/Pt sc." };
02503
02504 double thisMini[] = { 0.8, -0.1,
02505 0., 0., 0., 0.,
02506 1.0,
02507 -0.1, 0. };
02508 if( muonType == 1 ) {
02509 double thisMaxi[] = { 1., 1., 1., 1.,
02510 1., 1., 1., 1.,1.};
02511
02512 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02513 } else {
02514 double thisMaxi[] = { 2., 0.1,
02515 0.01, 0.1, 0.1, 1.,
02516 4.,
02517 0.1, 0.1 };
02518
02519 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02520 }
02521 }
02522 };
02523
02527
02528 template <class T>
02529 class resolutionFunctionType30 : public resolutionFunctionBase<T> {
02530 public:
02531 resolutionFunctionType30() { this->parNum_ = 27; }
02532
02533 virtual double sigmaPt(const double & pt, const double & eta, const T & parval)
02534 {
02535 double fabsEta = std::fabs(eta);
02536
02537 double ptPart = parval[13]*pt;
02538 if( fabsEta > 2.0 ) {
02539 ptPart = parval[22]*pt + parval[23]*pt*pt;
02540 }
02541 else if( fabsEta > 1.4 ) {
02542 ptPart = parval[20]*pt + parval[21]*pt*pt;
02543 }
02544 if(fabsEta<parval[0]) {
02545 return( ptPart + parval[1] + parval[2]*fabsEta + parval[3]*eta*eta );
02546 }
02547
02548 else if( fabsEta < parval[14] ) {
02549 double x_1 = parval[0];
02550 double y_1 = parval[1] + parval[2]*parval[0] + parval[3]*parval[0]*parval[0];
02551 double x_2 = parval[14];
02552 double y_2 = parval[4] + parval[5]*std::fabs((parval[14]-parval[7])) + parval[6]*(parval[14]-parval[7])*(parval[14]-parval[7]);
02553 return( (fabsEta - x_1)*(y_2 - y_1)/(x_2 - x_1) + y_1 );
02554 }
02555 else if( fabsEta < parval[15] ) {
02556 return( ptPart + parval[4] + parval[5]*std::fabs(fabsEta-parval[7]) + parval[6]*(fabsEta-parval[7])*(fabsEta-parval[7]) );
02557 }
02558 else {
02559 return( ptPart + parval[16] + parval[17]*std::fabs(fabsEta-parval[19]) + parval[18]*(fabsEta-parval[19])*(fabsEta-parval[19]) );
02560 }
02561 }
02562
02563 virtual double sigmaPtError(const double & pt, const double & eta, const T & parval, const T & parError)
02564 {
02565 double fabsEta = std::fabs(eta);
02566
02567 double ptPart = parError[13]*pt;
02568 if( fabsEta > 2.0 ) {
02569 ptPart = parError[22]*pt + parError[23]*pt*pt;
02570 }
02571 else if( fabsEta > 1.4 ) {
02572 ptPart = parError[20]*pt + parError[21]*pt*pt;
02573 }
02574 if(fabsEta<parval[0]) {
02575 return( ptPart + parError[1] + parError[2]*fabsEta + parError[3]*eta*eta );
02576 }
02577
02578 else if( fabsEta < parval[14] ) {
02579 double x_1 = parval[0];
02580 double y_1 = parval[1] + parval[2]*parval[0] + parval[3]*parval[0]*parval[0];
02581 double x_2 = parval[14];
02582 double y_2 = parval[4] + parval[5]*std::fabs((parval[14]-parval[7])) + parval[6]*(parval[14]-parval[7])*(parval[14]-parval[7]);
02583 double lineValue = (fabsEta - x_1)*(y_2 - y_1)/(x_2 - x_1) + y_1;
02584
02585
02586 y_1 = parval[1] + parError[1] + (parval[2] + parError[2])*parval[0] + (parval[3] + parError[3])*parval[0]*parval[0];
02587
02588 y_2 = parval[4] + parError[4] + (parval[5] + parError[5])*std::fabs((parval[14]-parval[7])) + (parval[6] + parError[6])*(parval[14]-parval[7])*(parval[14]-parval[7]);
02589 double lineValuePlusError = (fabsEta - x_1)*(y_2 - y_1)/(x_2 - x_1) + y_1;
02590
02591 return(lineValuePlusError - lineValue );
02592 }
02593 else if( fabsEta < parval[15] ) {
02594 return( ptPart + parError[4] + parError[5]*std::fabs(fabsEta-parval[7]) + parError[6]*(fabsEta-parval[7])*(fabsEta-parval[7]) );
02595 }
02596 else {
02597 return( ptPart + parError[16] + parError[17]*std::fabs(fabsEta-parval[19]) + parError[18]*(fabsEta-parval[19])*(fabsEta-parval[19]) );
02598 }
02599 }
02600
02601
02602 virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02603
02604 double fabsEta = std::fabs(eta);
02605 double value = parval[8] + parval[9]*fabsEta + parval[10]*eta*eta + parval[11]*fabsEta*fabsEta*fabsEta;
02606 if( value > 0 ) {
02607 return( value );
02608 }
02609 return 0;
02610 }
02611
02612
02613 virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02614 return( parval[12] );
02615 }
02616
02617 virtual double covPt1Pt2(const double & pt1, const double & eta1, const double & pt2, const double & eta2, const T & parval)
02618 {
02619 return parval[24] + std::fabs(pt1 - pt2)*parval[25] + std::fabs(eta1 - eta2)*parval[26];
02620 }
02621
02622 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const std::vector<int> & parResolOrder, const int muonType)
02623 {
02624 std::vector<ParameterSet> parSet(this->parNum_);
02625
02626 parSet[0] = ParameterSet( "etaTransition", 0.0001, 0., 2. );
02627 parSet[1] = ParameterSet( "constantCentral", 0.00001, 0., 1. );
02628 parSet[2] = ParameterSet( "linearEtaCentral", 0.00001, 0., 1. );
02629 parSet[3] = ParameterSet( "quadraticEtaCentral", 0.000001, 0., 1. );
02630 parSet[4] = ParameterSet( "constantForward", 0.00001, 0., 1. );
02631 parSet[5] = ParameterSet( "linearEtaForward", 0.00001, 0., 1. );
02632 parSet[6] = ParameterSet( "quadraticEtaForward", 0.000001, 0., 1. );
02633 parSet[7] = ParameterSet( "vertexForward", 0.0001, 0., 3. );
02634 parSet[8] = ParameterSet( "cotgThetaConstant", 0.00001, -1., 1. );
02635 parSet[9] = ParameterSet( "cotgThetaFactor", 0.00001, -1., 1. );
02636 parSet[10] = ParameterSet( "cotgThetaDenominatorTerm", 0.000001, -1., 1. );
02637 parSet[11] = ParameterSet( "cotgThetaLinearPt", 0.000001, -1., 1. );
02638 parSet[12] = ParameterSet( "sigmaPhi", 0.0001, 0., 1. );
02639 parSet[13] = ParameterSet( "barrelLinearPt", 0.00001, 0., 1. );
02640 parSet[14] = ParameterSet( "split", 0.0001, 0., 3. );
02641 parSet[15] = ParameterSet( "veryForwardSplit", 0.0001, 0., 3. );
02642 parSet[16] = ParameterSet( "constantVeryForward", 0.00001, 0., 1. );
02643 parSet[17] = ParameterSet( "linearEtaVeryForward", 0.00001, 0., 1. );
02644 parSet[18] = ParameterSet( "quadraticEtaVeryForward", 0.000001, 0., 1. );
02645 parSet[19] = ParameterSet( "vertexVeryForward", 0.0001, 0., 3. );
02646 parSet[20] = ParameterSet( "endcapsLinearPt", 0.00001, -1., 1. );
02647 parSet[21] = ParameterSet( "endcapsQuadraticPt", 0.000001, -1., 1. );
02648 parSet[22] = ParameterSet( "veryForwardLinearPt", 0.00001, -1., 1. );
02649 parSet[23] = ParameterSet( "veryForwardQuadraticPt", 0.000001, -1., 1. );
02650 parSet[24] = ParameterSet( "covPt1Pt2Constant", 0.000001, -1., 1. );
02651 parSet[25] = ParameterSet( "covPt1Pt2DeltaPt", 0.000001, -1., 1. );
02652 parSet[26] = ParameterSet( "covPt1Pt2DeltaEta", 0.000001, -1., 1. );
02653
02654 std::cout << "setting parameters" << std::endl;
02655 this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parSet );
02656 }
02657 };
02658
02659
02660
02661
02662
02663
02664
02666 resolutionFunctionBase<double *> * resolutionFunctionService( const int identifier );
02667
02669 resolutionFunctionBase<std::vector<double> > * resolutionFunctionVecService( const int identifier );
02670
02695 class backgroundFunctionBase {
02696 public:
02697 backgroundFunctionBase(const double & lowerLimit, const double & upperLimit) :
02698 lowerLimit_(lowerLimit), upperLimit_(upperLimit) {}
02699 virtual ~backgroundFunctionBase()
02700 {
02701 delete functionForIntegral_;
02702 };
02703 virtual double operator()( const double * parval, const double & mass, const double & eta ) const = 0;
02704 virtual int parNum() const { return parNum_; }
02706 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const std::vector<double>::const_iterator & parBgrIt, const std::vector<int>::const_iterator & parBgrOrderIt, const int muonType) = 0;
02707 virtual TF1* functionForIntegral(const std::vector<double>::const_iterator & parBgrIt) const
02708 {
02709 functionForIntegral_ = new FunctionForIntegral(this, parBgrIt);
02710 TF1 * backgroundFunctionForIntegral = new TF1("backgroundFunctionForIntegral", functionForIntegral_,
02711 lowerLimit_, upperLimit_, this->parNum_);
02712 return( backgroundFunctionForIntegral );
02713 }
02714 virtual double fracVsEta(const double * parval, const double & resEta) const { return 1.; }
02715
02716 protected:
02717 int parNum_;
02718 double lowerLimit_;
02719 double upperLimit_;
02721 virtual void setPar(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname,
02722 const std::vector<double>::const_iterator & parBgrIt, const std::vector<int>::const_iterator & parBgrOrderIt,
02723 double* thisStep, double* thisMini, double* thisMaxi, TString* thisParName ) {
02724 for( int iPar=0; iPar<this->parNum_; ++iPar ) {
02725 Start[iPar] = *(parBgrIt+iPar);
02726 Step[iPar] = thisStep[iPar];
02727 Mini[iPar] = thisMini[iPar];
02728 Maxi[iPar] = thisMaxi[iPar];
02729 ind[iPar] = *(parBgrOrderIt+iPar);
02730 parname[iPar] = thisParName[iPar];
02731 }
02732 }
02733 class FunctionForIntegral
02734 {
02735 public:
02736 FunctionForIntegral( const backgroundFunctionBase * function,
02737 const std::vector<double>::const_iterator & parBgrIt ) :
02738 function_(function)
02739 {
02740 parval_ = new double[function_->parNum()];
02741 for( int i=0; i < function_->parNum(); ++i ) {
02742 parval_[i] = *(parBgrIt+i);
02743 }
02744 }
02745 ~FunctionForIntegral()
02746 {
02747 delete parval_;
02748 }
02749 double operator()(const double * mass, const double *) const
02750 {
02751
02752 return( (*function_)(parval_, *mass, 0.) );
02753 }
02754 protected:
02755 const backgroundFunctionBase * function_;
02756 double * parval_;
02757 };
02758 mutable FunctionForIntegral * functionForIntegral_;
02759 };
02760
02762
02763 class backgroundFunctionType1 : public backgroundFunctionBase
02764 {
02765 public:
02772 backgroundFunctionType1(const double & lowerLimit, const double & upperLimit) :
02773 backgroundFunctionBase(lowerLimit, upperLimit)
02774 { this->parNum_ = 2; }
02775 virtual double operator()( const double * parval, const double & mass, const double & eta ) const
02776 {
02777 double a = 1.;
02778 double b = parval[1];
02779
02780 double norm = -(a*lowerLimit_ + b*lowerLimit_*lowerLimit_/2.);
02781
02782 if( -a/b > upperLimit_ ) norm += a*upperLimit_ + b*upperLimit_*upperLimit_/2.;
02783 else norm += -a*a/(2*b);
02784
02785 if( mass < -a/b && norm != 0 ) return (a + b*mass)/norm;
02786 else return 0;
02787 }
02788 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const std::vector<double>::const_iterator & parBgrIt, const std::vector<int>::const_iterator & parBgrOrderIt, const int muonType) {
02789 double thisStep[] = {0.01, 0.01};
02790 TString thisParName[] = {"Bgr fraction", "Constant", "Linear"};
02791 if( muonType == 1 ) {
02792 double thisMini[] = {0.0, -300.};
02793 double thisMaxi[] = {1.0, 0.};
02794 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02795 } else {
02796 double thisMini[] = {0.0, -300.};
02797 double thisMaxi[] = {1.0, 0.};
02798 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02799 }
02800 }
02801 };
02803
02804 class backgroundFunctionType2 : public backgroundFunctionBase {
02805 public:
02811 backgroundFunctionType2(const double & lowerLimit, const double & upperLimit) :
02812 backgroundFunctionBase(lowerLimit, upperLimit)
02813 { this->parNum_ = 2; }
02814 virtual double operator()( const double * parval, const double & mass, const double & eta ) const
02815 {
02816 double Bgrp2 = parval[1];
02817 double norm = -(exp(-Bgrp2*upperLimit_) - exp(-Bgrp2*lowerLimit_))/Bgrp2;
02818 if( norm != 0 ) return exp(-Bgrp2*mass)/norm;
02819 else return 0.;
02820 }
02821 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const std::vector<double>::const_iterator & parBgrIt, const std::vector<int>::const_iterator & parBgrOrderIt, const int muonType) {
02822 double thisStep[] = {0.01, 0.01};
02823 TString thisParName[] = {"Bgr fraction", "Bgr slope"};
02824 if( muonType == 1 ) {
02825 double thisMini[] = {0.0, 0.};
02826 double thisMaxi[] = {1.0, 10.};
02827 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02828 } else {
02829 double thisMini[] = {0.0, 0.};
02830 double thisMaxi[] = {1.0, 10.};
02831 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02832 }
02833 }
02834
02835
02836
02837
02838
02839
02840
02841
02842
02843
02844
02845
02846
02847 };
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885
02886
02887
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899
02901
02902 class backgroundFunctionType4 : public backgroundFunctionBase
02903 {
02904 public:
02910 backgroundFunctionType4(const double & lowerLimit, const double & upperLimit) :
02911 backgroundFunctionBase(lowerLimit, upperLimit)
02912 { this->parNum_ = 4; }
02913 virtual double operator()( const double * parval, const double & mass, const double & eta ) const
02914 {
02915 double Bgrp2 = parval[1] + parval[2]*eta*eta;
02916 double norm = -(exp(-Bgrp2*upperLimit_) - exp(-Bgrp2*lowerLimit_))/Bgrp2;
02917 if( norm != 0 ) return exp(-Bgrp2*mass)/norm;
02918 else return 0.;
02919 }
02920 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const std::vector<double>::const_iterator & parBgrIt, const std::vector<int>::const_iterator & parBgrOrderIt, const int muonType) {
02921 double thisStep[] = {0.01, 0.01, 0.01, 0.01};
02922 TString thisParName[] = {"Bgr fraction", "Bgr slope", "Bgr slope eta^2 dependence", "background fraction eta dependence"};
02923 if( muonType == 1 ) {
02924 double thisMini[] = {0.0, 0., 0., -1.};
02925 double thisMaxi[] = {1.0, 10., 10., 1.};
02926 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02927 } else {
02928 double thisMini[] = {0.0, 0., -1., -1.};
02929 double thisMaxi[] = {1.0, 10., 1., 1.};
02930 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02931 }
02932 }
02933 virtual double fracVsEta(const double * parval, const double & resEta) const
02934 {
02935 return( 1. - parval[3]*resEta*resEta );
02936 }
02937 };
02938
02940
02941 class backgroundFunctionType5 : public backgroundFunctionBase
02942 {
02943 public:
02948 backgroundFunctionType5(const double & lowerLimit, const double & upperLimit) :
02949 backgroundFunctionBase(lowerLimit, upperLimit)
02950 { this->parNum_ = 3; }
02951 virtual double operator()( const double * parval, const double & mass, const double & eta ) const
02952 {
02953 double b = parval[1];
02954
02955 double a = 1 + parval[2]*eta*eta;
02956
02957 double norm = -(a*lowerLimit_ + b*lowerLimit_*lowerLimit_/2.);
02958
02959 if( -a/b > upperLimit_ ) norm += a*upperLimit_ + b*upperLimit_*upperLimit_/2.;
02960 else norm += -a*a/(2*b);
02961
02962 if( mass < -a/b && norm != 0 ) return (a + b*mass)/norm;
02963 else return 0;
02964 }
02965 virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const std::vector<double>::const_iterator & parBgrIt, const std::vector<int>::const_iterator & parBgrOrderIt, const int muonType) {
02966 double thisStep[] = {0.01, 0.01, 0.01};
02967 TString thisParName[] = {"Bgr fraction", "Constant", "Linear"};
02968 if( muonType == 1 ) {
02969 double thisMini[] = {0.0, 0., -300.};
02970 double thisMaxi[] = {1.0, 300., 0.};
02971 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02972 } else {
02973 double thisMini[] = {0.0, 0., -300.};
02974 double thisMaxi[] = {1.0, 300., 0.};
02975 this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02976 }
02977 }
02978 };
02979
02980
02982 backgroundFunctionBase * backgroundFunctionService( const int identifier, const double & lowerLimit, const double & upperLimit );
02983
02984 #endif // FUNCTIONS_H