Heppy/interface/MuScleFitCorrector_Functions.h

0001 #ifndef PhysicsTools_Heppy_MuScleFitCorrector_Functions_h
0002 #define PhysicsTools_Heppy_MuScleFitCorrector_Functions_h
0003
0004 /**
0005  * Scale function classes
0006  * Author M. De Mattia - 18/11/2008
0007  * Author S. Casasso   - 25/10/2012
0008  * Author E. Migliore  - 25/10/2012
0009  */
0010
0011 #include <iostream>
0012 #include <vector>
0013 #include <cmath>
0014 #include "TMath.h"
0015 #include "TString.h"
0016 #include "TF1.h"
0017 #include "TRandom.h"
0018
0019 /**
0020  * Used to define parameters inside the functions.
0021  */
0022
0023 namespace heppy {
0024   struct ParSet {
0025     ParSet() {}
0026     ParSet(const TString& inputName, const double& inputStep, const double& inputMini, const double& inputMaxi)
0027         : step(inputStep), mini(inputMini), maxi(inputMaxi) {
0028       std::cout << "setting name = " << inputName << std::endl;
0029       name = inputName;
0030     }
0031     TString name;
0032     double step, mini, maxi;
0033   };
0034
0035   // ----------------------- //
0036   //     Scale functors      //
0037   // ----------------------- //
0038
0039   template <class T>
0040   class scaleFunctionBase {
0041   public:
0042     virtual double scale(
0043         const double& pt, const double& eta, const double& phi, const int chg, const T& parScale) const = 0;
0044     virtual ~scaleFunctionBase() = 0;
0045     virtual int parNum() const { return parNum_; }
0046
0047   protected:
0048     int parNum_;
0049     virtual void setPar(double* Start,
0050                         double* Step,
0051                         double* Mini,
0052                         double* Maxi,
0053                         int* ind,
0054                         TString* parname,
0055                         const T& parResol,
0056                         const std::vector<int>& parResolOrder,
0057                         const std::vector<ParSet>& parSet) {
0058       if (int(parSet.size()) != this->parNum_) {
0059         std::cout << "Error: wrong number of parameter initializations = " << parSet.size()
0060                   << ". Number of parameters is " << this->parNum_ << std::endl;
0061         exit(1);
0062       }
0063       for (int iPar = 0; iPar < this->parNum_; ++iPar) {
0064         Start[iPar] = parResol[iPar];
0065         Step[iPar] = parSet[iPar].step;
0066         Mini[iPar] = parSet[iPar].mini;
0067         Maxi[iPar] = parSet[iPar].maxi;
0068         ind[iPar] = parResolOrder[iPar];
0069         parname[iPar] = parSet[iPar].name;
0070       }
0071     }
0072   };
0073
0074   template <class T>
0075   inline scaleFunctionBase<T>::~scaleFunctionBase() {
0076   }  // defined even though it's pure virtual; should be faster this way.
0077
0078   //
0079   // Curvature: (linear eta + sinusoidal in phi (both in 5 eta bins)) * global scale
0080   // ------------------------------------------------------------
0081   template <class T>
0082   class scaleFunction50 : public scaleFunctionBase<T> {
0083   public:
0084     scaleFunction50() { this->parNum_ = 27; }
0085     virtual double scale(
0086         const double& pt, const double& eta, const double& phi, const int chg, const T& parScale) const {
0087       double ampl(0), phase(0), twist(0), ampl2(0), freq2(0), phase2(0);
0088
0089       // very bwd bin
0090       if (eta < parScale[4]) {
0091         ampl = parScale[1];
0092         phase = parScale[2];
0093         ampl2 = parScale[21];
0094         freq2 = parScale[22];
0095         phase2 = parScale[23];
0096         twist =
0097             parScale[3] * (eta - parScale[4]) + parScale[7] * (parScale[4] - parScale[8]) + parScale[11] * parScale[8];
0098         // bwd bin
0099       } else if (parScale[4] <= eta && eta < parScale[8]) {
0100         ampl = parScale[5];
0101         phase = parScale[6];
0102         twist = parScale[7] * (eta - parScale[8]) + parScale[11] * parScale[8];
0103         // barrel bin
0104       } else if (parScale[8] <= eta && eta < parScale[12]) {
0105         ampl = parScale[9];
0106         phase = parScale[10];
0107         twist = parScale[11] * eta;
0108         // fwd bin
0109       } else if (parScale[12] <= eta && eta < parScale[16]) {
0110         ampl = parScale[13];
0111         phase = parScale[14];
0112         twist = parScale[15] * (eta - parScale[12]) + parScale[11] * parScale[12];
0113         // very fwd bin
0114       } else if (parScale[16] < eta) {
0115         ampl = parScale[17];
0116         phase = parScale[18];
0117         ampl2 = parScale[24];
0118         freq2 = parScale[25];
0119         phase2 = parScale[26];
0120         twist = parScale[19] * (eta - parScale[16]) + parScale[15] * (parScale[16] - parScale[12]) +
0121                 parScale[11] * parScale[12];
0122       }
0123
0124       // apply the correction
0125       double curv = (1. + parScale[0]) * ((double)chg / pt - twist - ampl * sin(phi + phase) -
0126                                           ampl2 * sin(freq2 * phi + phase2) - 0.5 * parScale[20]);
0127       return 1. / ((double)chg * curv);
0128     }
0129   };
0130
0131   // ----------------------- //
0132   //   Resolution functors   //
0133   // ----------------------- //
0134
0135   template <class T>
0136   class resolutionFunctionBase {
0137   public:
0138     virtual double sigmaPt(const double& pt, const double& eta, const T& parval) = 0;
0139
0140     resolutionFunctionBase() {}
0141     virtual ~resolutionFunctionBase() = 0;
0142     virtual int parNum() const { return parNum_; }
0143
0144   protected:
0145     int parNum_;
0146   };
0147   template <class T>
0148   inline resolutionFunctionBase<T>::~resolutionFunctionBase() {
0149   }  // defined even though it's pure virtual; should be faster this way.
0150
0151   template <class T>
0152   class resolutionFunction45 : public resolutionFunctionBase<T> {
0153   public:
0154     resolutionFunction45() { this->parNum_ = 13; }
0155
0156     inline double getGEO(const double& pt, const double& eta, const T& parval) { return parval[0]; }
0157
0158     inline double getMS(const double& pt, const double& eta, const T& parval) {
0159       if (eta < -2.0)
0160         return (parval[1]);
0161       if (eta < -1.8)
0162         return (parval[2]);
0163       if (eta < -1.6)
0164         return (parval[3]);
0165       if (eta < -1.2)
0166         return (parval[4]);
0167       if (eta < -0.8)
0168         return (parval[5]);
0169       if (eta < 0.)
0170         return (parval[6]);
0171       if (eta < 0.8)
0172         return (parval[7]);
0173       if (eta < 1.2)
0174         return (parval[8]);
0175       if (eta < 1.6)
0176         return (parval[9]);
0177       if (eta < 1.8)
0178         return (parval[10]);
0179       if (eta < 2.0)
0180         return (parval[11]);
0181       return (parval[12]);
0182     }
0183
0184     virtual double sigmaPt(const double& pt, const double& eta, const T& parval) {
0185       return pt * getGEO(pt, eta, parval) + getMS(pt, eta, parval);
0186     }
0187   };
0188
0189   template <class T>
0190   class resolutionFunction46 : public resolutionFunctionBase<T> {
0191   public:
0192     resolutionFunction46() { this->parNum_ = 13; }
0193
0194     int etaBin(const double& eta) {
0195       // std::cout << "for eta = " << eta << ", bin = " << bin << std::endl;
0196
0197       if (eta < -2.0)
0198         return 1;
0199       if (eta < -1.8)
0200         return 2;
0201       if (eta < -1.6)
0202         return 3;
0203       if (eta < -1.2)
0204         return 4;
0205       if (eta < -0.8)
0206         return 5;
0207       if (eta < 0.)
0208         return 6;
0209       if (eta < 0.8)
0210         return 7;
0211       if (eta < 1.2)
0212         return 8;
0213       if (eta < 1.6)
0214         return 9;
0215       if (eta < 1.8)
0216         return 10;
0217       if (eta < 2.0)
0218         return 11;
0219       return 12;
0220     }
0221
0222     virtual double sigmaPt(const double& pt, const double& eta, const T& parval) {
0223       return sqrt(pow(parval[0] * pt, 2) + pow(parval[etaBin(eta)], 2));
0224     }
0225   };
0226
0227   // parametrization as sum in quadrature
0228   // Geometric and MSC both as function of eta, adding straight lines between parabolas wrt type51
0229   template <class T>
0230   class resolutionFunction57 : public resolutionFunctionBase<T> {
0231   public:
0232     resolutionFunction57() { this->parNum_ = 17; }
0233
0234     inline double getGEO(const double& pt, const double& eta, const T& parval) {
0235       // geometrical term
0236       double qGEO(0);
0237       if (eta < parval[0]) {
0238         qGEO = parval[10] * (eta - parval[0]) * (eta - parval[0]) + parval[11];
0239       } else if (parval[0] <= eta && eta < parval[3]) {
0240         qGEO = parval[11];
0241       } else {
0242         qGEO = parval[12] * (eta - parval[3]) * (eta - parval[3]) + parval[11];
0243       }
0244       return qGEO;
0245     }
0246
0247     inline double centralParabola(const double& pt, const double& eta, const T& parval) {
0248       return parval[4] + parval[5] * eta * eta;
0249     }
0250
0251     inline double middleParabola(const double& pt, const double& eta, const T& parval) {
0252       return parval[15] + parval[16] * eta * eta;
0253     }
0254
0255     inline double leftParabola(const double& pt, const double& eta, const T& parval) {
0256       return parval[6] + parval[7] * (eta - parval[0]) * (eta - parval[0]);
0257     }
0258
0259     inline double rightParabola(const double& pt, const double& eta, const T& parval) {
0260       return parval[8] + parval[9] * (eta - parval[3]) * (eta - parval[3]);
0261     }
0262
0263     inline double leftLine(const double& pt, const double& eta, const T& parval) {
0264       double x_1 = parval[13];
0265       double y_1 = middleParabola(pt, parval[13], parval);
0266       double x_2 = parval[0];
0267       double y_2 = leftParabola(pt, parval[0], parval);
0268       return ((eta - x_1) * (y_2 - y_1) / (x_2 - x_1) + y_1);
0269     }
0270
0271     inline double rightLine(const double& pt, const double& eta, const T& parval) {
0272       double x_1 = parval[14];
0273       double y_1 = middleParabola(pt, parval[14], parval);
0274       double x_2 = parval[3];
0275       double y_2 = rightParabola(pt, parval[3], parval);
0276       return ((eta - x_1) * (y_2 - y_1) / (x_2 - x_1) + y_1);
0277     }
0278
0279     inline double getMSC(const double& pt, const double& eta, const T& parval) {
0280       // MSC term
0281       double qMSC(0);
0282       if (eta < parval[0]) {
0283         qMSC = leftParabola(pt, eta, parval);
0284       } else if (parval[0] <= eta && eta < parval[13]) {
0285         qMSC = leftLine(pt, eta, parval);
0286       } else if (parval[13] <= eta && eta < parval[1]) {
0287         qMSC = middleParabola(pt, eta, parval);
0288       } else if (parval[1] <= eta && eta < parval[2]) {
0289         qMSC = centralParabola(pt, eta, parval);
0290       } else if (parval[2] <= eta && eta < parval[14]) {
0291         qMSC = middleParabola(pt, eta, parval);
0292       } else if (parval[14] <= eta && eta < parval[3]) {
0293         qMSC = rightLine(pt, eta, parval);
0294       } else {
0295         qMSC = rightParabola(pt, eta, parval);
0296       }
0297       return qMSC;
0298     }
0299
0300     virtual double sigmaPt(const double& pt, const double& eta, const T& parval) {
0301       double qGEO = getGEO(pt, eta, parval);
0302       double qMSC = getMSC(pt, eta, parval);
0303       return sqrt(pow(qGEO * pt, 2) + pow(qMSC, 2));
0304     }
0305   };
0306
0307   // Service to build the scale functor corresponding to the passed identifier
0308   scaleFunctionBase<double*>* scaleFunctionService(const int identifier) {
0309     switch (identifier) {
0310       case (50):
0311         return (new scaleFunction50<double*>);
0312         break;
0313       default:
0314         std::cout << "scaleFunctionService error: wrong identifier = " << identifier << std::endl;
0315         exit(1);
0316     }
0317   }
0318
0319   // Service to build the resolution functor corresponding to the passed identifier
0320   resolutionFunctionBase<double*>* resolutionFunctionService(const int identifier) {
0321     switch (identifier) {
0322       case (45):
0323         return (new resolutionFunction45<double*>);
0324         break;
0325       case (46):
0326         return (new resolutionFunction46<double*>);
0327         break;
0328       case (57):
0329         return (new resolutionFunction57<double*>);
0330         break;
0331       default:
0332         std::cout << "resolutionFunctService error: wrong identifier = " << identifier << std::endl;
0333         exit(1);
0334     }
0335   }
0336
0337 }  // namespace heppy
0338
0339 #endif