CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h

Go to the documentation of this file.
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 // Bias and scale functors //
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() { }  // defined even though it's pure virtual; should be faster this way.
00074 // No scale
00075 // --------
00076 template <class T>
00077 class scaleFunctionType0 : public scaleFunctionBase<T> {
00078 public:
00079   scaleFunctionType0() {
00080     // One of the two is required. This follows from when templates are used by the compiler and the names lookup rules in c++.
00081     // scaleFunctionBase<T>::parNum_ = 0;
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 // Linear in pt
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   // Fill the scaleVec with neutral parameters
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 // Linear in |eta|
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   // Fill the scaleVec with neutral parameters
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 // Sinusoidal in phi
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt and |eta|
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt and sinusoidal in phi
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   // Fill the scaleVec with neutral parameters
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 // Linear in |eta| and sinusoidal in phi
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt and |eta| and sinusoidal in phi
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt and parabolic in |eta|
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   // Fill the scaleVec with neutral parameters
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       // double thisMini[] = {0.985, -0.002, -0.005, -0.005};
00332       // double thisMaxi[] = {1.015, 0.002, 0.005, 0.005};
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 // Exponential in pt
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   // Fill the scaleVec with neutral parameters
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 // Parabolic in pt
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt, sinusoidal in phi with muon sign
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt, parabolic in eta, sinusoidal in phi with muon sign
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt, parabolic in eta, sinusoidal in phi with muon sign
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     // else {
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   // Fill the scaleVec with neutral parameters
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 // linear in pt and cubic in |eta|
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 //     for( int i=0; i<10; ++i ) {
00484 //       std::cout << " parScale["<<i<<"] = " << parScale[i];
00485 //     }
00486 //     std::cout << "   newPt = " << ( parScale[0] +
00487 //                                parScale[1]*pt + parScale[2]*pt*pt + parScale[3]*pt*pt*pt +
00488 //                                parScale[4]*std::fabs(eta) + parScale[5]*eta*eta + parScale[6]*std::fabs(eta*eta*eta) +
00489 //                                parScale[7]*eta*eta*eta*eta + parScale[8]*std::fabs(eta*eta*eta*eta*eta) +
00490 //                                parScale[9]*eta*eta*eta*eta*eta*eta )*pt << std::endl;
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   // Fill the scaleVec with neutral parameters
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 // linear in pt and cubic in |eta|
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   // Fill the scaleVec with neutral parameters
00546   virtual void resetParameters(std::vector<double> * scaleVec) const {
00547     // For the first, reset to the original pt region separator
00548     scaleVec->push_back(originalPtRegionSeparator_);
00549     // The next two are the scale in the two regions
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   // Fill the scaleVec with neutral parameters
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   // Fill the scaleVec with neutral parameters
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   // Fill the scaleVec with neutral parameters
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 // ---- R.C.Nov.09 ---
00697 // Scale function for Z mass (misalignment STARTUP scenario) corrections
00698 // Linear in pt, sinusoidal in phi (muon-charge dependent) and parabolic in Eta
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   // Fill the scaleVec with neutral parameters
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 // This function allows to use three different pt functions for three pt ranges
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   // Fill the scaleVec with neutral parameters
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 // Linear in pt and up to cubic in |eta| with possible eta asymmetry: two parabolic branches are used one for eta+ and one for eta-
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   // Fill the scaleVec with neutral parameters
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 // Function built to correct STARTUP MC
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     // Set to 0: use the same parameters for negative and positive muons
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   // Fill the scaleVec with neutral parameters
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                          // , 0.0001, 0.01, 0.0001, 0.01};
00877     TString thisParName[] = {"Phi offset",
00878                              "amplitude pos phi", "phase pos phi",
00879                              "amplitude neg phi", "phase neg phi"};
00880                              // "amplitude pos charge pos phi", "phase pos charge pos phi",
00881                              // "amplitude pos charge neg phi", "phase pos charge neg phi",
00882                              // "amplitude neg charge pos phi", "phase neg charge pos phi",
00883                              // "amplitude neg charge neg phi", "phase neg charge neg phi" };
00884     if( muonType == 1 ) {
00885       double thisMini[] = {0.9, -0.3, -0.3, -0.3, -0.3};
00886                            // , -0.3, -0.3, -0.3, -0.3};
00887       double thisMaxi[] = {1.1,  0.3,  0.3,  0.3,  0.3};
00888                            // , 0.3,  0.3,  0.3,  0.3};
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                            // , -0.1, -3, -0.1, -3};
00893       double thisMaxi[] = {1.1,   0.1,  3,  0.1,  3};
00894                            // , 0.1,  3,  0.1,  3};
00895       this->setPar( Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName );
00896     }
00897   }
00898 };
00899 
00900 
00901 // Function built to correct STARTUP MC
00902 // Independent parameters for mu+ and mu-
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     // Set to 0: use the same parameters for negative and positive muons
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   // Fill the scaleVec with neutral parameters
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                              // "amplitude pos phi", "phase pos phi",
00968                              // "amplitude neg phi", "phase neg phi"};
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 // Function built to correct STARTUP MC
01004 // Independent parameters for mu+ and mu-
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     // Set to 0: use the same parameters for negative and positive muons
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   // Fill the scaleVec with neutral parameters
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                              // "amplitude pos phi", "phase pos phi",
01082                              // "amplitude neg phi", "phase neg phi"};
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 // Function built to correct STARTUP MC
01114 // Analytical description in eta
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     // Set to 0: use the same parameters for negative and positive muons
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   // Fill the scaleVec with neutral parameters
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                              // "amplitude pos phi", "phase pos phi",
01175                              // "amplitude neg phi", "phase neg phi"};
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 // Built for the first 100/nb of J/Psi in data
01216 // It has eta dependent corrections only for |eta| > parScale[6] and separate parabolic corrections for eta > 0 or < 0. 
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   // Fill the scaleVec with neutral parameters
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 // Built for the first 100/nb of J/Psi in data
01264 // It has eta dependent corrections only for |eta| > parScale[6] and separate parabolic corrections for eta > 0 or < 0. 
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   // Fill the scaleVec with neutral parameters
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 // Smear functors //
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() { }  // defined even though it's pure virtual; should be faster this way.
01368 
01369 // No smearing
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 // The 3 parameters of smearType1 are: pt dependence of pt smear, phi smear and
01376 // cotgtheta smear.
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 // The six parameters of SmearType=4 are respectively:
01410 // Pt dep. of Pt res., |eta| dep. of Pt res., Phi res., |eta| res.,
01411 // |eta| dep. of |eta| res., Pt^2 dep. of Pt res.
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 //Smearing for MC correction based on the resolution function Type 15 for misaligned MC
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){//eta in right endcap
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{//eta in left endcap
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     // ATTENTION: This point has a big error and it is very displaced from the rest of the distribution.
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 // // Defined globally...
01508 // static smearFunctionBase * smearFunctionArray[] = {
01509 //   new smearFunctionType0,
01510 //   new smearFunctionType1,
01511 //   new smearFunctionType2,
01512 //   new smearFunctionType3,
01513 //   new smearFunctionType4,
01514 //   new smearFunctionType5
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() { }  // defined even though it's pure virtual; should be faster this way.
01572 
01573 // Resolution Type 1
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 // Resolution Type 6
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 // Resolution Type 7
01640 template <class T>
01641 class resolutionFunctionType7 : public resolutionFunctionBase<T> {
01642  public:
01643   resolutionFunctionType7() { this->parNum_ = 12; }
01644   // linear in pt and quadratic in eta
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   // 1/pt in pt and quadratic in eta
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   // 1/pt in pt and quadratic in eta
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 // Resolution Type 8
01681 template <class T>
01682 class resolutionFunctionType8 : public resolutionFunctionBase<T> {
01683  public:
01684   resolutionFunctionType8() { this->parNum_ = 12; }
01685   // linear in pt and by points in eta
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   // 1/pt in pt and quadratic in eta
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   // 1/pt in pt and quadratic in eta
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 // Resolution Type 9
01745 template <class T>
01746 class resolutionFunctionType9 : public resolutionFunctionBase<T> {
01747  public:
01748   resolutionFunctionType9() { this->parNum_ = 31; }
01749   // linear in pt and by points in eta
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   // 1/pt in pt and quadratic in eta
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   // 1/pt in pt and quadratic in eta
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     // More detailed taken from Upsilon2S
01837 //     if( 0.0 < eta && eta <= 0.1 ) return( (0.006496598 + 0.006713836)/2 );
01838 //     else if( 0.1 < eta && eta <= 0.2 ) return( (0.006724315 + 0.006787474)/2 );
01839 //     else if( 0.2 < eta && eta <= 0.3 ) return( (0.007284029 + 0.007293643)/2 );
01840 //     else if( 0.3 < eta && eta <= 0.4 ) return( (0.008138282 + 0.008187387)/2 );
01841 //     else if( 0.4 < eta && eta <= 0.5 ) return( (0.008174111 + 0.008030496)/2 );
01842 //     else if( 0.5 < eta && eta <= 0.6 ) return( (0.008126558 + 0.008100443)/2 );
01843 //     else if( 0.6 < eta && eta <= 0.7 ) return( (0.008602069 + 0.008626195)/2 );
01844 //     else if( 0.7 < eta && eta <= 0.8 ) return( (0.009187699 + 0.009090244)/2 );
01845 //     else if( 0.8 < eta && eta <= 0.9 ) return( (0.009835283 + 0.009875661)/2 );
01846 //     else if( 0.9 < eta && eta <= 1.0 ) return( (0.01156847 + 0.011774)/2);
01847 //     else if( 1.0 < eta && eta <= 1.1 ) return( (0.01319311 + 0.01312528)/2 );
01848 //     else if( 1.1 < eta && eta <= 1.2 ) return( (0.01392963 + 0.01413793)/2 );
01849 //     else if( 1.2 < eta && eta <= 1.3 ) return( (0.01430238 + 0.01385749)/2 );
01850 //     else if( 1.3 < eta && eta <= 1.4 ) return( (0.01409375 + 0.01450355)/2 );
01851 //     else if( 1.4 < eta && eta <= 1.5 ) return( (0.01395235 + 0.01419122)/2 );
01852 //     else if( 1.5 < eta && eta <= 1.6 ) return( (0.01384032 + 0.01354162)/2 );
01853 //     else if( 1.6 < eta && eta <= 1.7 ) return( (0.01325593 + 0.01302663)/2 );
01854 //     else if( 1.7 < eta && eta <= 1.8 ) return( (0.01365382 + 0.01361993)/2 );
01855 //     else if( 1.8 < eta && eta <= 1.9 ) return( (0.01516075 + 0.01514115)/2 );
01856 //     else if( 1.9 < eta && eta <= 2.0 ) return( (0.01587837 + 0.01561742)/2 );
01857 //     else if( 2.0 < eta && eta <= 2.1 ) return( (0.01696865 + 0.01760318)/2 );
01858 //     else if( 2.1 < eta && eta <= 2.2 ) return( (0.01835451 + 0.01887852)/2 );
01859 //     else if( 2.2 < eta && eta <= 2.3 ) return( (0.02116863 + 0.02254953)/2 );
01860 //     else if( 2.3 < eta && eta <= 2.4 ) return( (0.0224906 + 0.02158211)/2 );
01861 
01862     // Less detailed
01863 //     if( 0.0 < eta && eta <= 0.2 ) return( (0.006496598 + 0.006713836 + 0.006724315 + 0.006787474)/4 );
01864 //     else if( 0.2 < eta && eta <= 0.4 ) return( (0.007284029 + 0.007293643 + 0.008138282 + 0.008187387)/4 );
01865 //     else if( 0.4 < eta && eta <= 0.6 ) return( (0.008174111 + 0.008030496 + 0.008126558 + 0.008100443)/4 );
01866 //     else if( 0.6 < eta && eta <= 0.8 ) return( (0.008602069 + 0.008626195 + 0.009187699 + 0.009090244)/4 );
01867 //     else if( 0.8 < eta && eta <= 1.0 ) return( (0.009835283 + 0.009875661 + 0.01156847 + 0.011774)/4 );
01868 //     else if( 1.0 < eta && eta <= 1.2 ) return( (0.01319311 + 0.01312528 + 0.01392963 + 0.01413793)/4 );
01869 //     else if( 1.2 < eta && eta <= 1.4 ) return( (0.01430238 + 0.01385749 + 0.01409375 + 0.01450355)/4 );
01870 //     else if( 1.4 < eta && eta <= 1.6 ) return( (0.01395235 + 0.01419122 + 0.01384032 + 0.01354162)/4 );
01871 //     else if( 1.6 < eta && eta <= 1.8 ) return( (0.01325593 + 0.01302663 + 0.01365382 + 0.01361993)/4 );
01872 //     else if( 1.8 < eta && eta <= 2.0 ) return( (0.01516075 + 0.01514115 + 0.01587837 + 0.01561742)/4 );
01873 //     else if( 2.0 < eta && eta <= 2.2 ) return( (0.01696865 + 0.01760318 + 0.01835451 + 0.01887852)/4 );
01874 //     // else if( 2.2 < eta && eta <= 2.4 ) return( (0.02116863 + 0.02254953 + 0.0224906 + 0.02158211)/4 );
01875 
01876 //     return ( border );
01877 
01878     // From MuonGun
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 // Resolution Type 10
01899 template <class T>
01900 class resolutionFunctionType10 : public resolutionFunctionBase<T> {
01901  public:
01902   resolutionFunctionType10() { this->parNum_ = 21; }
01903   // linear in pt and by points in eta
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   // 1/pt in pt and quadratic in eta
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   // 1/pt in pt and quadratic in eta
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 // Resolution Type 11
01956 template <class T>
01957 class resolutionFunctionType11 : public resolutionFunctionBase<T> {
01958  public:
01959   resolutionFunctionType11() { this->parNum_ = 8; }
01960   // linear in pt and by points in eta
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   // 1/pt in pt and quadratic in eta
01969   virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
01970     return( 0.004 );
01971   }
01972   // 1/pt in pt and quadratic in eta
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 // Resolution Type 12
02001 template <class T>
02002 class resolutionFunctionType12 : public resolutionFunctionBase<T> {
02003  public:
02004   resolutionFunctionType12() { this->parNum_ = 15; }
02005   // linear in pt and by points in eta
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       // To impose continuity we require that the parval[0] of type11 is
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   // 1/pt in pt and quadratic in eta
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   // // 1/pt in pt and quadratic in eta
02027   // virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02028   //   return( parval[15]+parval[16]/pt + parval[17]*std::fabs(eta)+parval[18]*eta*eta );
02029   // }
02030 
02031   // constant sigmaCotgTh
02032   // virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02033   //   return( 0.004 );
02034   // }
02035 
02036   // constant sigmaPhi
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                           // 0.00002, 0.0002, 0.00000002, 0.000002 };
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                               // "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc." };
02053     double thisMini[] = { 0.8, -1.1, 0., -1.1,
02054                           0., 0.0005, 0.0005, 0.001,
02055                           1.4, 0., 0.,
02056                           // -0.001, 0.002, -0.0001, -0.0001 };
02057                           -0.1, 0., -0.1, -0.1 };
02058                           // -0.0001, 0.0005, -0.0001, -0.00001 };
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                             // 1., 1., 1., 1. };
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                             // 0.001, 0.005, 0.00004, 0.0007 };
02072                             0.1, 0.05, 0.1, 0.1 };
02073                             // 0.001, 0.01, -0.0000015, 0.0004 };
02074       this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName );
02075     }
02076   }
02077 };
02078 
02084 // Resolution Type 13
02085 template <class T>
02086 class resolutionFunctionType13 : public resolutionFunctionBase<T> {
02087  public:
02088   resolutionFunctionType13() { this->parNum_ = 15; }
02089   // linear in pt and by points in eta
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       // To impose continuity we require that
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         // return( parval[15] + ptPart + parval[4]*fabsEta + parval[5]*eta*eta );
02102       }
02103       // return( par2+ ptPart + parval[10]*std::fabs((fabsEta-parval[12])) + parval[11]*(fabsEta-parval[12])*(fabsEta-parval[12]) );
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   // 1/pt in pt and quadratic in eta
02109   virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02110     return( 0.004 );
02111   }
02112   // 1/pt in pt and quadratic in eta
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 //                          0.001 };
02123     TString thisParName[] = { "etaTransition1", "offsetEtaHigh", "coeffOverPt", "coeffHighPt",
02124                               "linearEtaCentral", "parabEtaCentral", "linearEtaHighEta", "parabEtaHighEta",
02125                               "centerEtaHighEta", "etaTransition2", "linearEtaSecondEta", "parabEtaSecondEta",
02126                               "centerEtaSecondEta", "linearPt", "quadraticPt" };
02127 //                              "scale2" };
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 //                          -1.0 };
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 //                            1. };
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 //                            2.0 };
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   // linear in pt and by points in eta
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   // 1/pt in pt and quadratic in eta
02165   virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02166     return( 0.004 );
02167   }
02168   // 1/pt in pt and quadratic in eta
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 // par6,par7,par10
02190 // Resolution Type 15. For misaligned data. Linear in pt, parabolic in eta, regions separated: barrl+overlap, right endcap, left endcap.
02191 template <class T>
02192 class resolutionFunctionType15 : public resolutionFunctionBase<T> {
02193  public:
02194   resolutionFunctionType15() { this->parNum_ = 7; }
02195   // linear in pt and parabolic in eta
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) {//eta in barrel + overlap
02203       double par = - parval[1]*0.6*0.6;
02204       return( par + ptPart + parval[1]*eta*eta );
02205     }
02206     else if (eta>1.3){//eta in right endcap
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{//eta in left endcap
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   // 1/pt in pt and quadratic in eta
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   // 1/pt in pt and quadratic in eta
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 //Divides the resolution vs pt in three regions, barrel, overlap and endcaps. It gives excellent results for at least 15K Z
02255 template <class T>
02256 class resolutionFunctionType17 : public resolutionFunctionBase<T> {
02257  public:
02258   resolutionFunctionType17() { this->parNum_ = 18; }
02259   // linear in pt and parabolic in eta
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) {//eta in barrel
02266       return( ptPartBar + parval[3] + parval[4]*fabsEta);
02267     }
02268     else if( (eta > 0.9 && eta <= 1.4) || (eta < -0.9 && eta > -1.4)){ //eta in overlap
02269       return( ptPartOvlap + parval[3] + parval[4]*eta + parval[5]*eta*eta);
02270     }
02271     else if (eta>1.4){//eta in right endcap
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 {//eta in left endcap
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   // 1/pt in pt and quadratic in eta
02281   virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02282     return( 0.004 );
02283   }
02284   // 1/pt in pt and quadratic in eta
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   // linear in pt and parabolic in eta
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)) {//eta in barrel + overlap
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){//eta in right endcap
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{//eta in left endcap
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   // 1/pt in pt and quadratic in eta
02353   virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02354     return( 0.004 );
02355   }
02356   // 1/pt in pt and quadratic in eta
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 // Resolution Type 19
02395 // Same as type 8, but the sigmaPhi and sigmaCotgTh are not free. This way the function results as having less parameters.
02396 // This was done to verify if fixed parameters have an influence in the computation of errors by minuit.
02397 template <class T>
02398 class resolutionFunctionType19 : public resolutionFunctionBase<T> {
02399  public:
02400   resolutionFunctionType19() { this->parNum_ = 4; }
02401   // linear in pt and by points in eta
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   // 1/pt in pt and quadratic in eta
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   // 1/pt in pt and quadratic in eta
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   // linear in pt and by points in eta
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       // To impose continuity we require that the parval[0] of type11 is
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   // 1/pt in pt and quadratic in eta
02473   virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02474     return( parval[7]+parval[8]/pt  );
02475   }
02476 
02477   // // 1/pt in pt and quadratic in eta
02478   // virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) {
02479   //   return( parval[15]+parval[16]/pt + parval[17]*std::fabs(eta)+parval[18]*eta*eta );
02480   // }
02481 
02482   // constant sigmaCotgTh
02483   // virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02484   //   return( 0.004 );
02485   // }
02486 
02487   // constant sigmaPhi
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                           // 0.00002, 0.0002, 0.00000002, 0.000002 };
02499     TString thisParName[] = { "etaTransition", "offsetEtaHigh", 
02500                               "linaerEtaCentral", "parabEtaCentral", "linaerEtaHigh", "parabEtaHigh",
02501                               "secondParabolaCenter",
02502                               "Cth res. sc.", "Cth res. 1/Pt sc." };
02503                               // "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc." };
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 // Resolution Type 30
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     // Return a line connecting the two parabolas
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     // Return a line connecting the two parabolas
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       // x_1 = parval[0];
02586       y_1 = parval[1] + parError[1] + (parval[2] + parError[2])*parval[0] + (parval[3] + parError[3])*parval[0]*parval[0];
02587       // x_2 = parval[14];
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   // 1/pt in pt and quadratic in eta
02602   virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) {
02603     // return( parval[8] + parval[9]/(pt + parval[10]) + parval[11]*pt );
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   // constant sigmaPhi
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     // name, step, mini, maxi
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 // ------------ ATTENTION ----------- //
02662 // Other functions are not in for now //
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       // FIXME: this is a gross approximation. The function should be integrated in eta over the sample.
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   // virtual double fracVsEta(const double * parval, const double & resEta) const
02839   // {
02840   //   // return( 0.6120 - 0.0225*eta*eta );
02841   //   return( 1. - 0.0225*resEta*resEta ); // so that a = 1 for eta = 0.
02842   // }
02843 
02844 
02845 
02846 
02847 };
02852 //
02853 //class backgroundFunctionType3 : public backgroundFunctionBase {
02854 // public:
02855 //  // pass parval[shift]
02856 //  backgroundFunctionType3(const double & lowerLimit, const double & upperLimit) :
02857 //    backgroundFunctionBase(lowerLimit, upperLimit)
02858 //    { this->parNum_ = 3; }
02859 //  virtual double operator()( const double * parval, const int resTotNum, const int ires, const bool * resConsidered,
02860 //                             const double * ResMass, const double ResHalfWidth[], const int MuonType, const double & mass, const int nbins ) {
02861 //    double PB = 0.;
02862 //    double Bgrp2 = parval[1];
02863 //    double Bgrp3 = parval[2];
02864 //    for (int ires=0; ires<resTotNum; ires++) {
02865 //      // In this case, by integrating between A and B, we get for f=exp(a-bx)+k:
02866 //      // INT = exp(a)/b*(exp(-bA)-exp(-bB))+k*(B-A) so our function, which in 1000 bins between A and B
02867 //      // gets a total of 1, is f = (exp(a-bx)+k)*(B-A)/nbins / (INT)
02868 //      // ----------------------------------------------------------------------------------------------
02869 //      if (resConsidered[ires]) {
02870 //      if (exp(-Bgrp2*(ResMass[ires]-ResHalfWidth[ires]))-exp(-Bgrp2*(ResMass[ires]+ResHalfWidth[ires]))>0) {
02871 //        PB += (exp(-Bgrp2*mass)+Bgrp3) *
02872 //          2*ResHalfWidth[ires]/(double)nbins /
02873 //          ( (exp(-Bgrp2*(ResMass[ires]-ResHalfWidth[ires]))-exp(-Bgrp2*(ResMass[ires]+ResHalfWidth[ires])))/
02874 //            Bgrp2 + Bgrp3*2*ResHalfWidth[ires] );
02875 //      } else {
02876 //        std::cout << "Impossible to compute Background probability! - some fix needed - Bgrp2=" << Bgrp2 << std::endl;
02877 //      }
02878 //      }
02879 //    }
02880 //    return PB;
02881 //  }
02882 //  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) {
02883 //    double thisStep[] = {0.1, 0.001, 0.1};
02884 //    TString thisParName[] = {"Bgr fraction", "Bgr slope", "Bgr constant"};
02885 //    if( muonType == 1 ) {
02886 //      double thisMini[] = {0.0, 0.000000001, 0.0};
02887 //      double thisMaxi[] = {1.0, 0.2, 1000};
02888 //      this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02889 //    } else {
02890 //      double thisMini[] = {0.0, 0.000000001, 0.0};
02891 //      double thisMaxi[] = {1.0, 0.2, 1000};
02892 //      this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
02893 //    }
02894 //  }
02895 //  virtual TF1* functionForIntegral(const std::vector<double>::const_iterator & parBgrIt) const {return 0;};
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 ); // so that a = 1 for eta = 0.
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     // double c = parval[2];
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