Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:31:30

0001 #ifndef FUNCTIONS_H
0002 #define FUNCTIONS_H
0003 
0004 #include <iostream>
0005 #include <vector>
0006 #include <cmath>
0007 #include "TMath.h"
0008 #include "TString.h"
0009 #include "TF1.h"
0010 #include "TRandom.h"
0011 #include "MuonAnalysis/MomentumScaleCalibration/interface/SigmaPtDiff.h"
0012 
0013 /**
0014  * Used to define parameters inside the functions.
0015  */
0016 struct ParameterSet {
0017   ParameterSet() {}
0018   ParameterSet(const TString& inputName, const double& inputStep, const double& inputMini, const double& inputMaxi)
0019       : step(inputStep), mini(inputMini), maxi(inputMaxi) {
0020     std::cout << "setting name = " << inputName << std::endl;
0021     name = inputName;
0022   }
0023   TString name;
0024   double step, mini, maxi;
0025 };
0026 
0027 // ----------------------- //
0028 // Bias and scale functors //
0029 // ----------------------- //
0030 /** The correct functor is selected at job start in the constructor.
0031  * The pt value is taken by reference and modified internally.
0032  * eta, phi and chg are taken by const reference.<br>
0033  * Made into a template so that it can be used with arrays too
0034  * (parval for the scale fit is an array, because Lykelihood is an
0035  * extern C function, because TMinuit asks it).<br>
0036  * Note that in the array case it takes the pointer by const reference,
0037  * thus the elements of the array are modifiable.
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   /// This method is used to reset the scale parameters to neutral values (useful for iterations > 0)
0046   virtual void resetParameters(std::vector<double>* scaleVec) const {
0047     std::cout << "ERROR: the resetParameters method must be defined in each scale function" << std::endl;
0048     std::cout << "Please add it to the scaleFunction you are using" << std::endl;
0049     exit(1);
0050   }
0051   /// This method is used to differentiate parameters among the different functions
0052   virtual void setParameters(double* Start,
0053                              double* Step,
0054                              double* Mini,
0055                              double* Maxi,
0056                              int* ind,
0057                              TString* parname,
0058                              const T& parScale,
0059                              const std::vector<int>& parScaleOrder,
0060                              const int muonType) = 0;
0061   virtual void setParameters(double* Start,
0062                              double* Step,
0063                              double* Mini,
0064                              double* Maxi,
0065                              int* ind,
0066                              TString* parname,
0067                              const T& parResol,
0068                              const std::vector<int>& parResolOrder,
0069                              const std::vector<double>& parStep,
0070                              const std::vector<double>& parMin,
0071                              const std::vector<double>& parMax,
0072                              const int muonType) {
0073     std::cout << "The method setParameters must be implemented for this scale function" << std::endl;
0074     exit(1);
0075   }
0076   virtual int parNum() const { return parNum_; }
0077 
0078 protected:
0079   int parNum_;
0080   /// This method sets the parameters
0081   virtual void setPar(double* Start,
0082                       double* Step,
0083                       double* Mini,
0084                       double* Maxi,
0085                       int* ind,
0086                       TString* parname,
0087                       const T& parScale,
0088                       const std::vector<int>& parScaleOrder,
0089                       double* thisStep,
0090                       double* thisMini,
0091                       double* thisMaxi,
0092                       TString* thisParName) {
0093     for (int iPar = 0; iPar < this->parNum_; ++iPar) {
0094       Start[iPar] = parScale[iPar];
0095       Step[iPar] = thisStep[iPar];
0096       Mini[iPar] = thisMini[iPar];
0097       Maxi[iPar] = thisMaxi[iPar];
0098       ind[iPar] = parScaleOrder[iPar];
0099       parname[iPar] = thisParName[iPar];
0100     }
0101   }
0102   virtual void setPar(double* Start,
0103                       double* Step,
0104                       double* Mini,
0105                       double* Maxi,
0106                       int* ind,
0107                       TString* parname,
0108                       const T& parResol,
0109                       const std::vector<int>& parResolOrder,
0110                       const std::vector<ParameterSet>& parSet) {
0111     if (int(parSet.size()) != this->parNum_) {
0112       std::cout << "Error: wrong number of parameter initializations = " << parSet.size()
0113                 << ". Number of parameters is " << this->parNum_ << std::endl;
0114       exit(1);
0115     }
0116     for (int iPar = 0; iPar < this->parNum_; ++iPar) {
0117       Start[iPar] = parResol[iPar];
0118       Step[iPar] = parSet[iPar].step;
0119       Mini[iPar] = parSet[iPar].mini;
0120       Maxi[iPar] = parSet[iPar].maxi;
0121       ind[iPar] = parResolOrder[iPar];
0122       parname[iPar] = parSet[iPar].name;
0123     }
0124   }
0125 };
0126 
0127 // SLIMMED VERSION
0128 // The full set of scaleFunction developed in the past can be found at
0129 // https://github.com/scasasso/cmssw/blob/test_binned_function/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h#L100-L5177
0130 
0131 template <class T>
0132 inline scaleFunctionBase<T>::~scaleFunctionBase() {
0133 }  // defined even though it's pure virtual; should be faster this way.
0134 
0135 // No scale
0136 // --------
0137 template <class T>
0138 class scaleFunctionType0 : public scaleFunctionBase<T> {
0139 public:
0140   scaleFunctionType0() {
0141     // One of the two is required. This follows from when templates are used by the compiler and the names lookup rules in c++.
0142     // scaleFunctionBase<T>::parNum_ = 0;
0143     this->parNum_ = 0;
0144   }
0145   double scale(const double& pt, const double& eta, const double& phi, const int chg, const T& parScale) const override {
0146     return pt;
0147   }
0148   void resetParameters(std::vector<double>* scaleVec) const override {}
0149   void setParameters(double* Start,
0150                      double* Step,
0151                      double* Mini,
0152                      double* Maxi,
0153                      int* ind,
0154                      TString* parname,
0155                      const T& parScale,
0156                      const std::vector<int>& parScaleOrder,
0157                      const int muonType) override {}
0158 };
0159 
0160 //
0161 // Curvature: (linear eta + sinusoidal in phi (both in 5 eta bins)) * global scale
0162 // ------------------------------------------------------------
0163 template <class T>
0164 class scaleFunctionType50 : public scaleFunctionBase<T> {
0165 public:
0166   scaleFunctionType50() { this->parNum_ = 27; }
0167   double scale(const double& pt, const double& eta, const double& phi, const int chg, const T& parScale) const override {
0168     double ampl(0), phase(0), twist(0), ampl2(0), freq2(0), phase2(0);
0169 
0170     // very bwd bin
0171     if (eta < parScale[4]) {
0172       ampl = parScale[1];
0173       phase = parScale[2];
0174       ampl2 = parScale[21];
0175       freq2 = parScale[22];
0176       phase2 = parScale[23];
0177       twist =
0178           parScale[3] * (eta - parScale[4]) + parScale[7] * (parScale[4] - parScale[8]) + parScale[11] * parScale[8];
0179       // bwd bin
0180     } else if (parScale[4] <= eta && eta < parScale[8]) {
0181       ampl = parScale[5];
0182       phase = parScale[6];
0183       twist = parScale[7] * (eta - parScale[8]) + parScale[11] * parScale[8];
0184       // barrel bin
0185     } else if (parScale[8] <= eta && eta < parScale[12]) {
0186       ampl = parScale[9];
0187       phase = parScale[10];
0188       twist = parScale[11] * eta;
0189       // fwd bin
0190     } else if (parScale[12] <= eta && eta < parScale[16]) {
0191       ampl = parScale[13];
0192       phase = parScale[14];
0193       twist = parScale[15] * (eta - parScale[12]) + parScale[11] * parScale[12];
0194       // very fwd bin
0195     } else if (parScale[16] < eta) {
0196       ampl = parScale[17];
0197       phase = parScale[18];
0198       ampl2 = parScale[24];
0199       freq2 = parScale[25];
0200       phase2 = parScale[26];
0201       twist = parScale[19] * (eta - parScale[16]) + parScale[15] * (parScale[16] - parScale[12]) +
0202               parScale[11] * parScale[12];
0203     }
0204 
0205     // apply the correction
0206     double curv = (1. + parScale[0]) * ((double)chg / pt - twist - ampl * sin(phi + phase) -
0207                                         ampl2 * sin((int)freq2 * phi + phase2) - 0.5 * parScale[20]);
0208     return 1. / ((double)chg * curv);
0209   }
0210   // Fill the scaleVec with neutral parameters
0211   void resetParameters(std::vector<double>* scaleVec) const override {
0212     //    scaleVec->push_back(1);
0213     for (int i = 0; i < this->parNum_; ++i) {
0214       scaleVec->push_back(0);
0215     }
0216   }
0217   void setParameters(double* Start,
0218                      double* Step,
0219                      double* Mini,
0220                      double* Maxi,
0221                      int* ind,
0222                      TString* parname,
0223                      const T& parScale,
0224                      const std::vector<int>& parScaleOrder,
0225                      const int muonType) override {
0226     double thisStep[] = {0.000001, 0.000001, 0.01,     0.000001, 0.01,     0.000001, 0.01,     0.000001, 0.01,
0227                          0.000001, 0.01,     0.000001, 0.01,     0.000001, 0.01,     0.000001, 0.01,     0.000001,
0228                          0.01,     0.000001, 0.000001, 0.000001, 0.01,     0.01,     0.000001, 0.01,     0.01};
0229 
0230     TString thisParName[] = {"Curv global scale",
0231                              "Phi ampl eta vbwd",
0232                              "Phi phase eta vbwd",
0233                              "Twist eta vbwd",
0234                              "vbwd/bwd bndry",
0235                              "Phi ampl eta  bwd",
0236                              "Phi phase eta  bwd",
0237                              "Twist eta  bwd",
0238                              "bwd/bar bndry",
0239                              "Phi ampl eta  bar",
0240                              "Phi phase eta  bar",
0241                              "Twist eta  bar",
0242                              "bar/fwd bndry",
0243                              "Phi ampl eta  fwd",
0244                              "Phi phase eta  fwd",
0245                              "Twist eta  fwd",
0246                              "fwd/vfwd bndry",
0247                              "Phi ampl eta vfwd",
0248                              "Phi phase eta vfwd",
0249                              "Twist eta vfwd",
0250                              "Charge depend bias",
0251                              "Phi ampl eta vbwd (2nd harmon.)",
0252                              "Phi freq. eta vbwd (2nd harmon.)",
0253                              "Phi phase eta vbwd (2nd harmon.)",
0254                              "Phi ampl eta vfwd (2nd harmon.)",
0255                              "Phi freq. eta vfwd (2nd harmon.)",
0256                              "Phi phase eta vfwd (2nd harmon.)"};
0257 
0258     if (muonType == 1) {
0259       double thisMini[] = {-0.1,    -0.3,    -3.1416, -0.3, -2.6, -0.3,    -3.1416, -0.3, -2.6,
0260                            -0.3,    -3.1416, -0.3,    0.,   -0.3, -3.1416, -0.3,    0.,   -0.3,
0261                            -3.1416, -0.3,    -0.1,    -0.3, 1.,   -3.1416, -0.3,    1.,   -3.1416};
0262       double thisMaxi[] = {0.1,    0.3, 3.1416, 0.3, 0.,     0.3, 3.1416, 0.3, 0., 0.3,    3.1416, 0.3, 2.6,   0.3,
0263                            3.1416, 0.3, 2.6,    0.3, 3.1416, 0.3, 0.1,    0.3, 5., 3.1416, 0.3,    5.,  3.1416};
0264       this->setPar(
0265           Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName);
0266     } else {
0267       double thisMini[] = {-0.1,    -0.3,    -3.1416, -0.3, -2.6, -0.3,    -3.1416, -0.3, -2.6,
0268                            -0.3,    -3.1416, -0.3,    0.,   -0.3, -3.1416, -0.3,    0.,   -0.3,
0269                            -3.1416, -0.3,    -0.1,    -0.3, 1.,   -3.1416, -0.3,    1.,   -3.1416};
0270       double thisMaxi[] = {0.1,    0.3, 3.1416, 0.3, 0.,     0.3, 3.1416, 0.3, 0., 0.3,    3.1416, 0.3, 2.6,   0.3,
0271                            3.1416, 0.3, 2.6,    0.3, 3.1416, 0.3, 0.1,    0.3, 5., 3.1416, 0.3,    5.,  3.1416};
0272       this->setPar(
0273           Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName);
0274     }
0275   }
0276   void setParameters(double* Start,
0277                      double* Step,
0278                      double* Mini,
0279                      double* Maxi,
0280                      int* ind,
0281                      TString* parname,
0282                      const T& parScale,
0283                      const std::vector<int>& parScaleOrder,
0284                      const std::vector<double>& parStep,
0285                      const std::vector<double>& parMin,
0286                      const std::vector<double>& parMax,
0287                      const int muonType) override {
0288     if ((int(parStep.size()) != this->parNum_) || (int(parMin.size()) != this->parNum_) ||
0289         (int(parMax.size()) != this->parNum_)) {
0290       std::cout << "Error: par step or min or max do not match with number of parameters" << std::endl;
0291       std::cout << "parNum = " << this->parNum_ << std::endl;
0292       std::cout << "parStep.size() = " << parStep.size() << std::endl;
0293       std::cout << "parMin.size() = " << parMin.size() << std::endl;
0294       std::cout << "parMax.size() = " << parMax.size() << std::endl;
0295       exit(1);
0296     }
0297     std::vector<ParameterSet> parSet(this->parNum_);
0298     // name, step, mini, maxi
0299     parSet[0] = ParameterSet("Curv global scale", parStep[0], parMin[0], parMax[0]);
0300     parSet[1] = ParameterSet("Phi ampl  vbwd", parStep[1], parMin[1], parMax[1]);
0301     parSet[2] = ParameterSet("Phi phase vbwd", parStep[2], parMin[2], parMax[2]);
0302     parSet[3] = ParameterSet("Twist vbwd", parStep[3], parMin[3], parMax[3]);
0303     parSet[4] = ParameterSet("vbwd/bwd bndry", parStep[4], parMin[4], parMax[4]);
0304     parSet[5] = ParameterSet("Phi ampl   bwd", parStep[5], parMin[5], parMax[5]);
0305     parSet[6] = ParameterSet("Phi phase  bwd", parStep[6], parMin[6], parMax[6]);
0306     parSet[7] = ParameterSet("Twist  bwd", parStep[7], parMin[7], parMax[7]);
0307     parSet[8] = ParameterSet("bwd/bar bndry", parStep[8], parMin[8], parMax[8]);
0308     parSet[9] = ParameterSet("Phi ampl   bar", parStep[9], parMin[9], parMax[9]);
0309     parSet[10] = ParameterSet("Phi phase  bar", parStep[10], parMin[10], parMax[10]);
0310     parSet[11] = ParameterSet("Twist  bar", parStep[11], parMin[11], parMax[11]);
0311     parSet[12] = ParameterSet("bar/fwd bndry", parStep[12], parMin[12], parMax[12]);
0312     parSet[13] = ParameterSet("Phi ampl   fwd", parStep[13], parMin[13], parMax[13]);
0313     parSet[14] = ParameterSet("Phi phase  fwd", parStep[14], parMin[14], parMax[14]);
0314     parSet[15] = ParameterSet("Twist  fwd", parStep[15], parMin[15], parMax[15]);
0315     parSet[16] = ParameterSet("fwd/vfwd bndry", parStep[16], parMin[16], parMax[16]);
0316     parSet[17] = ParameterSet("Phi ampl  vfwd", parStep[17], parMin[17], parMax[17]);
0317     parSet[18] = ParameterSet("Phi phase vfwd", parStep[18], parMin[18], parMax[18]);
0318     parSet[19] = ParameterSet("Twist vfwd", parStep[19], parMin[19], parMax[19]);
0319     parSet[20] = ParameterSet("Charge depend bias", parStep[20], parMin[20], parMax[20]);
0320     parSet[21] = ParameterSet("Phi ampl vbwd (2nd)", parStep[21], parMin[21], parMax[21]);
0321     parSet[22] = ParameterSet("Phi freq vbwd (2nd)", parStep[22], parMin[22], parMax[22]);
0322     parSet[23] = ParameterSet("Phi phase vbwd (2nd)", parStep[23], parMin[23], parMax[23]);
0323     parSet[24] = ParameterSet("Phi ampl vfwd (2nd)", parStep[24], parMin[24], parMax[24]);
0324     parSet[25] = ParameterSet("Phi freq vfwd (2nd)", parStep[25], parMin[25], parMax[25]);
0325     parSet[26] = ParameterSet("Phi phase vfwd (2nd)", parStep[26], parMin[26], parMax[26]);
0326 
0327     std::cout << "setting parameters" << std::endl;
0328     for (int i = 0; i < this->parNum_; ++i) {
0329       std::cout << "parStep[" << i << "] = " << parStep[i] << ", parMin[" << i << "] = " << parMin[i] << ", parMax["
0330                 << i << "] = " << parMin[i] << std::endl;
0331     }
0332     this->setPar(Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, parSet);
0333   }
0334 };
0335 
0336 //
0337 // Curvature: binned function w/ constraints on eta rings
0338 // good results on Legacy 2011 MC
0339 // phi bins = [8,8,8,8,8,8,8] in 7 eta bins
0340 // ------------------------------------------------------------
0341 template <class T>
0342 class scaleFunctionType64 : public scaleFunctionBase<T> {
0343 public:
0344   scaleFunctionType64() { this->parNum_ = 50; }
0345   double scale(const double& pt, const double& eta, const double& phi, const int chg, const T& parScale) const override {
0346     double deltaK(0);
0347     double p0 = parScale[0];
0348 
0349     if (eta < -2.1 && eta > -2.4 && phi < -2.35619449019 && phi > -3.14159265359)
0350       deltaK = parScale[1];
0351     else if (eta < -2.1 && eta > -2.4 && phi < -1.57079632679 && phi > -2.35619449019)
0352       deltaK = parScale[2];
0353     else if (eta < -2.1 && eta > -2.4 && phi < -0.785398163397 && phi > -1.57079632679)
0354       deltaK = parScale[3];
0355     else if (eta < -2.1 && eta > -2.4 && phi < 0.0 && phi > -0.785398163397)
0356       deltaK = parScale[4];
0357     else if (eta < -2.1 && eta > -2.4 && phi < 0.785398163397 && phi > 0.0)
0358       deltaK = parScale[5];
0359     else if (eta < -2.1 && eta > -2.4 && phi < 1.57079632679 && phi > 0.785398163397)
0360       deltaK = parScale[6];
0361     else if (eta < -2.1 && eta > -2.4 && phi < 2.35619449019 && phi > 1.57079632679)
0362       deltaK = parScale[7];
0363     else if (eta < -2.1 && eta > -2.4 && phi < 3.14159265359 && phi > 2.35619449019)
0364       deltaK = -(parScale[1] + parScale[2] + parScale[3] + parScale[4] + parScale[5] + parScale[6] + parScale[7]);
0365 
0366     else if (eta < -1.5 && eta > -2.1 && phi < -2.35619449019 && phi > -3.14159265359)
0367       deltaK = parScale[8];
0368     else if (eta < -1.5 && eta > -2.1 && phi < -1.57079632679 && phi > -2.35619449019)
0369       deltaK = parScale[9];
0370     else if (eta < -1.5 && eta > -2.1 && phi < -0.785398163397 && phi > -1.57079632679)
0371       deltaK = parScale[10];
0372     else if (eta < -1.5 && eta > -2.1 && phi < 0.0 && phi > -0.785398163397)
0373       deltaK = parScale[11];
0374     else if (eta < -1.5 && eta > -2.1 && phi < 0.785398163397 && phi > 0.0)
0375       deltaK = parScale[12];
0376     else if (eta < -1.5 && eta > -2.1 && phi < 1.57079632679 && phi > 0.785398163397)
0377       deltaK = parScale[13];
0378     else if (eta < -1.5 && eta > -2.1 && phi < 2.35619449019 && phi > 1.57079632679)
0379       deltaK = parScale[14];
0380     else if (eta < -1.5 && eta > -2.1 && phi < 3.14159265359 && phi > 2.35619449019)
0381       deltaK = -(parScale[8] + parScale[9] + parScale[10] + parScale[11] + parScale[12] + parScale[13] + parScale[14]);
0382 
0383     else if (eta < -0.9 && eta > -1.5 && phi < -2.35619449019 && phi > -3.14159265359)
0384       deltaK = parScale[15];
0385     else if (eta < -0.9 && eta > -1.5 && phi < -1.57079632679 && phi > -2.35619449019)
0386       deltaK = parScale[16];
0387     else if (eta < -0.9 && eta > -1.5 && phi < -0.785398163397 && phi > -1.57079632679)
0388       deltaK = parScale[17];
0389     else if (eta < -0.9 && eta > -1.5 && phi < 0.0 && phi > -0.785398163397)
0390       deltaK = parScale[18];
0391     else if (eta < -0.9 && eta > -1.5 && phi < 0.785398163397 && phi > 0.0)
0392       deltaK = parScale[19];
0393     else if (eta < -0.9 && eta > -1.5 && phi < 1.57079632679 && phi > 0.785398163397)
0394       deltaK = parScale[20];
0395     else if (eta < -0.9 && eta > -1.5 && phi < 2.35619449019 && phi > 1.57079632679)
0396       deltaK = parScale[21];
0397     else if (eta < -0.9 && eta > -1.5 && phi < 3.14159265359 && phi > 2.35619449019)
0398       deltaK =
0399           -(parScale[15] + parScale[16] + parScale[17] + parScale[18] + parScale[19] + parScale[20] + parScale[21]);
0400 
0401     else if (eta < 0.9 && eta > -0.9 && phi < -2.35619449019 && phi > -3.14159265359)
0402       deltaK = parScale[22];
0403     else if (eta < 0.9 && eta > -0.9 && phi < -1.57079632679 && phi > -2.35619449019)
0404       deltaK = parScale[23];
0405     else if (eta < 0.9 && eta > -0.9 && phi < -0.785398163397 && phi > -1.57079632679)
0406       deltaK = parScale[24];
0407     else if (eta < 0.9 && eta > -0.9 && phi < 0.0 && phi > -0.785398163397)
0408       deltaK = parScale[25];
0409     else if (eta < 0.9 && eta > -0.9 && phi < 0.785398163397 && phi > 0.0)
0410       deltaK = parScale[26];
0411     else if (eta < 0.9 && eta > -0.9 && phi < 1.57079632679 && phi > 0.785398163397)
0412       deltaK = parScale[27];
0413     else if (eta < 0.9 && eta > -0.9 && phi < 2.35619449019 && phi > 1.57079632679)
0414       deltaK = parScale[28];
0415     else if (eta < 0.9 && eta > -0.9 && phi < 3.14159265359 && phi > 2.35619449019)
0416       deltaK =
0417           -(parScale[22] + parScale[23] + parScale[24] + parScale[25] + parScale[26] + parScale[27] + parScale[28]);
0418 
0419     else if (eta < 1.5 && eta > 0.9 && phi < -2.35619449019 && phi > -3.14159265359)
0420       deltaK = parScale[29];
0421     else if (eta < 1.5 && eta > 0.9 && phi < -1.57079632679 && phi > -2.35619449019)
0422       deltaK = parScale[30];
0423     else if (eta < 1.5 && eta > 0.9 && phi < -0.785398163397 && phi > -1.57079632679)
0424       deltaK = parScale[31];
0425     else if (eta < 1.5 && eta > 0.9 && phi < 0.0 && phi > -0.785398163397)
0426       deltaK = parScale[32];
0427     else if (eta < 1.5 && eta > 0.9 && phi < 0.785398163397 && phi > 0.0)
0428       deltaK = parScale[33];
0429     else if (eta < 1.5 && eta > 0.9 && phi < 1.57079632679 && phi > 0.785398163397)
0430       deltaK = parScale[34];
0431     else if (eta < 1.5 && eta > 0.9 && phi < 2.35619449019 && phi > 1.57079632679)
0432       deltaK = parScale[35];
0433     else if (eta < 1.5 && eta > 0.9 && phi < 3.14159265359 && phi > 2.35619449019)
0434       deltaK =
0435           -(parScale[29] + parScale[30] + parScale[31] + parScale[32] + parScale[33] + parScale[34] + parScale[35]);
0436 
0437     else if (eta < 2.1 && eta > 1.5 && phi < -2.35619449019 && phi > -3.14159265359)
0438       deltaK = parScale[36];
0439     else if (eta < 2.1 && eta > 1.5 && phi < -1.57079632679 && phi > -2.35619449019)
0440       deltaK = parScale[37];
0441     else if (eta < 2.1 && eta > 1.5 && phi < -0.785398163397 && phi > -1.57079632679)
0442       deltaK = parScale[38];
0443     else if (eta < 2.1 && eta > 1.5 && phi < 0.0 && phi > -0.785398163397)
0444       deltaK = parScale[39];
0445     else if (eta < 2.1 && eta > 1.5 && phi < 0.785398163397 && phi > 0.0)
0446       deltaK = parScale[40];
0447     else if (eta < 2.1 && eta > 1.5 && phi < 1.57079632679 && phi > 0.785398163397)
0448       deltaK = parScale[41];
0449     else if (eta < 2.1 && eta > 1.5 && phi < 2.35619449019 && phi > 1.57079632679)
0450       deltaK = parScale[42];
0451     else if (eta < 2.1 && eta > 1.5 && phi < 3.14159265359 && phi > 2.35619449019)
0452       deltaK =
0453           -(parScale[36] + parScale[37] + parScale[38] + parScale[39] + parScale[40] + parScale[41] + parScale[42]);
0454 
0455     else if (eta < 2.4 && eta > 2.1 && phi < -2.35619449019 && phi > -3.14159265359)
0456       deltaK = parScale[43];
0457     else if (eta < 2.4 && eta > 2.1 && phi < -1.57079632679 && phi > -2.35619449019)
0458       deltaK = parScale[44];
0459     else if (eta < 2.4 && eta > 2.1 && phi < -0.785398163397 && phi > -1.57079632679)
0460       deltaK = parScale[45];
0461     else if (eta < 2.4 && eta > 2.1 && phi < 0.0 && phi > -0.785398163397)
0462       deltaK = parScale[46];
0463     else if (eta < 2.4 && eta > 2.1 && phi < 0.785398163397 && phi > 0.0)
0464       deltaK = parScale[47];
0465     else if (eta < 2.4 && eta > 2.1 && phi < 1.57079632679 && phi > 0.785398163397)
0466       deltaK = parScale[48];
0467     else if (eta < 2.4 && eta > 2.1 && phi < 2.35619449019 && phi > 1.57079632679)
0468       deltaK = parScale[49];
0469     else if (eta < 2.4 && eta > 2.1 && phi < 3.14159265359 && phi > 2.35619449019)
0470       deltaK =
0471           -(parScale[43] + parScale[44] + parScale[45] + parScale[46] + parScale[47] + parScale[48] + parScale[49]);
0472     else {
0473       std::cout << "This should really not happen, this muon has eta = " << eta << "and phi = " << phi << std::endl;
0474       exit(1);
0475     }
0476 
0477     // apply the correction
0478     double curv = (double)chg / pt;
0479     return 1. / ((double)chg * (1 + p0) * (curv + deltaK));
0480   }
0481   // Fill the scaleVec with neutral parameters
0482   void resetParameters(std::vector<double>* scaleVec) const override {
0483     //    scaleVec->push_back(1);
0484     for (int i = 0; i < this->parNum_; ++i) {
0485       scaleVec->push_back(0);
0486     }
0487   }
0488   void setParameters(double* Start,
0489                      double* Step,
0490                      double* Mini,
0491                      double* Maxi,
0492                      int* ind,
0493                      TString* parname,
0494                      const T& parScale,
0495                      const std::vector<int>& parScaleOrder,
0496                      const int muonType) override {
0497     double thisStep[] = {0.000001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,
0498                          0.001,    0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,
0499                          0.001,    0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,
0500                          0.001,    0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,
0501                          0.001,    0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001};
0502 
0503     TString thisParName[] = {"Curv global scale", "deltaK bin1",  "deltaK bin2",  "deltaK bin3",  "deltaK bin4",
0504                              "deltaK bin5",       "deltaK bin6",  "deltaK bin7",  "deltaK bin8",  "deltaK bin9",
0505                              "deltaK bin10",      "deltaK bin11", "deltaK bin12", "deltaK bin13", "deltaK bin14",
0506                              "deltaK bin15",      "deltaK bin16", "deltaK bin17", "deltaK bin18", "deltaK bin19",
0507                              "deltaK bin20",      "deltaK bin21", "deltaK bin22", "deltaK bin23", "deltaK bin24",
0508                              "deltaK bin25",      "deltaK bin26", "deltaK bin27", "deltaK bin28", "deltaK bin29",
0509                              "deltaK bin30",      "deltaK bin31", "deltaK bin32", "deltaK bin33", "deltaK bin34",
0510                              "deltaK bin35",      "deltaK bin36", "deltaK bin37", "deltaK bin38", "deltaK bin39",
0511                              "deltaK bin40",      "deltaK bin41", "deltaK bin42", "deltaK bin43", "deltaK bin44",
0512                              "deltaK bin45",      "deltaK bin46", "deltaK bin47", "deltaK bin48", "deltaK bin49"};
0513 
0514     if (muonType == 1) {
0515       double thisMini[] = {-0.1,   -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0516                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0517                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0518                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0519                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005};
0520       double thisMaxi[] = {0.1,   0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005,
0521                            0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005,
0522                            0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005,
0523                            0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005};
0524       this->setPar(
0525           Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName);
0526     } else {
0527       double thisMini[] = {-0.1,   -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0528                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0529                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0530                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005,
0531                            -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005, -0.005};
0532       double thisMaxi[] = {0.1,   0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005,
0533                            0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005,
0534                            0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005,
0535                            0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005};
0536       this->setPar(
0537           Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, thisStep, thisMini, thisMaxi, thisParName);
0538     }
0539   }
0540   void setParameters(double* Start,
0541                      double* Step,
0542                      double* Mini,
0543                      double* Maxi,
0544                      int* ind,
0545                      TString* parname,
0546                      const T& parScale,
0547                      const std::vector<int>& parScaleOrder,
0548                      const std::vector<double>& parStep,
0549                      const std::vector<double>& parMin,
0550                      const std::vector<double>& parMax,
0551                      const int muonType) override {
0552     if ((int(parStep.size()) != this->parNum_) || (int(parMin.size()) != this->parNum_) ||
0553         (int(parMax.size()) != this->parNum_)) {
0554       std::cout << "Error: par step or min or max do not match with number of parameters" << std::endl;
0555       std::cout << "parNum = " << this->parNum_ << std::endl;
0556       std::cout << "parStep.size() = " << parStep.size() << std::endl;
0557       std::cout << "parMin.size() = " << parMin.size() << std::endl;
0558       std::cout << "parMax.size() = " << parMax.size() << std::endl;
0559       exit(1);
0560     }
0561     std::vector<ParameterSet> parSet(this->parNum_);
0562     // name, step, mini, maxi
0563     parSet[0] = ParameterSet("Curv global scale", parStep[0], parMin[0], parMax[0]);
0564     parSet[1] = ParameterSet("deltaK bin1", parStep[1], parMin[1], parMax[1]);
0565     parSet[2] = ParameterSet("deltaK bin2", parStep[2], parMin[2], parMax[2]);
0566     parSet[3] = ParameterSet("deltaK bin3", parStep[3], parMin[3], parMax[3]);
0567     parSet[4] = ParameterSet("deltaK bin4", parStep[4], parMin[4], parMax[4]);
0568     parSet[5] = ParameterSet("deltaK bin5", parStep[5], parMin[5], parMax[5]);
0569     parSet[6] = ParameterSet("deltaK bin6", parStep[6], parMin[6], parMax[6]);
0570     parSet[7] = ParameterSet("deltaK bin7", parStep[7], parMin[7], parMax[7]);
0571     parSet[8] = ParameterSet("deltaK bin9", parStep[8], parMin[8], parMax[8]);
0572     parSet[9] = ParameterSet("deltaK bin10", parStep[9], parMin[9], parMax[9]);
0573     parSet[10] = ParameterSet("deltaK bin11", parStep[10], parMin[10], parMax[10]);
0574     parSet[11] = ParameterSet("deltaK bin12", parStep[11], parMin[11], parMax[11]);
0575     parSet[12] = ParameterSet("deltaK bin13", parStep[12], parMin[12], parMax[12]);
0576     parSet[13] = ParameterSet("deltaK bin14", parStep[13], parMin[13], parMax[13]);
0577     parSet[14] = ParameterSet("deltaK bin15", parStep[14], parMin[14], parMax[14]);
0578     parSet[15] = ParameterSet("deltaK bin17", parStep[15], parMin[15], parMax[15]);
0579     parSet[16] = ParameterSet("deltaK bin18", parStep[16], parMin[16], parMax[16]);
0580     parSet[17] = ParameterSet("deltaK bin19", parStep[17], parMin[17], parMax[17]);
0581     parSet[18] = ParameterSet("deltaK bin20", parStep[18], parMin[18], parMax[18]);
0582     parSet[19] = ParameterSet("deltaK bin21", parStep[19], parMin[19], parMax[19]);
0583     parSet[20] = ParameterSet("deltaK bin22", parStep[20], parMin[20], parMax[20]);
0584     parSet[21] = ParameterSet("deltaK bin23", parStep[21], parMin[21], parMax[21]);
0585     parSet[22] = ParameterSet("deltaK bin25", parStep[22], parMin[22], parMax[22]);
0586     parSet[23] = ParameterSet("deltaK bin26", parStep[23], parMin[23], parMax[23]);
0587     parSet[24] = ParameterSet("deltaK bin27", parStep[24], parMin[24], parMax[24]);
0588     parSet[25] = ParameterSet("deltaK bin28", parStep[25], parMin[25], parMax[25]);
0589     parSet[26] = ParameterSet("deltaK bin29", parStep[26], parMin[26], parMax[26]);
0590     parSet[27] = ParameterSet("deltaK bin30", parStep[27], parMin[27], parMax[27]);
0591     parSet[28] = ParameterSet("deltaK bin31", parStep[28], parMin[28], parMax[28]);
0592     parSet[29] = ParameterSet("deltaK bin33", parStep[29], parMin[29], parMax[29]);
0593     parSet[30] = ParameterSet("deltaK bin34", parStep[30], parMin[30], parMax[30]);
0594     parSet[31] = ParameterSet("deltaK bin35", parStep[31], parMin[31], parMax[31]);
0595     parSet[32] = ParameterSet("deltaK bin36", parStep[32], parMin[32], parMax[32]);
0596     parSet[33] = ParameterSet("deltaK bin37", parStep[33], parMin[33], parMax[33]);
0597     parSet[34] = ParameterSet("deltaK bin38", parStep[34], parMin[34], parMax[34]);
0598     parSet[35] = ParameterSet("deltaK bin39", parStep[35], parMin[35], parMax[35]);
0599     parSet[36] = ParameterSet("deltaK bin41", parStep[36], parMin[36], parMax[36]);
0600     parSet[37] = ParameterSet("deltaK bin42", parStep[37], parMin[37], parMax[37]);
0601     parSet[38] = ParameterSet("deltaK bin43", parStep[38], parMin[38], parMax[38]);
0602     parSet[39] = ParameterSet("deltaK bin44", parStep[39], parMin[39], parMax[39]);
0603     parSet[40] = ParameterSet("deltaK bin45", parStep[40], parMin[40], parMax[40]);
0604     parSet[41] = ParameterSet("deltaK bin46", parStep[41], parMin[41], parMax[41]);
0605     parSet[42] = ParameterSet("deltaK bin47", parStep[42], parMin[42], parMax[42]);
0606     parSet[43] = ParameterSet("deltaK bin49", parStep[43], parMin[43], parMax[43]);
0607     parSet[44] = ParameterSet("deltaK bin50", parStep[44], parMin[44], parMax[44]);
0608     parSet[45] = ParameterSet("deltaK bin51", parStep[45], parMin[45], parMax[45]);
0609     parSet[46] = ParameterSet("deltaK bin52", parStep[46], parMin[46], parMax[46]);
0610     parSet[47] = ParameterSet("deltaK bin53", parStep[47], parMin[47], parMax[47]);
0611     parSet[48] = ParameterSet("deltaK bin54", parStep[48], parMin[48], parMax[48]);
0612     parSet[49] = ParameterSet("deltaK bin55", parStep[49], parMin[49], parMax[49]);
0613 
0614     std::cout << "setting parameters" << std::endl;
0615     for (int i = 0; i < this->parNum_; ++i) {
0616       std::cout << "parStep[" << i << "] = " << parStep[i] << ", parMin[" << i << "] = " << parMin[i] << ", parMax["
0617                 << i << "] = " << parMin[i] << std::endl;
0618     }
0619     this->setPar(Start, Step, Mini, Maxi, ind, parname, parScale, parScaleOrder, parSet);
0620   }
0621 };
0622 
0623 /// Service to build the scale functor corresponding to the passed identifier
0624 scaleFunctionBase<double*>* scaleFunctionService(const int identifier);
0625 
0626 /// Service to build the scale functor corresponding to the passed identifier when receiving a std::vector<double>
0627 scaleFunctionBase<std::vector<double> >* scaleFunctionVecService(const int identifier);
0628 
0629 // -------------- //
0630 // Smear functors //
0631 // -------------- //
0632 
0633 class smearFunctionBase {
0634 public:
0635   virtual void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) = 0;
0636   smearFunctionBase() {
0637     cotgth_ = 0.;
0638     gRandom_ = new TRandom();
0639   }
0640   virtual ~smearFunctionBase() = 0;
0641 
0642 protected:
0643   void smearEta(double& eta) {
0644     double theta;
0645     if (cotgth_ != 0) {
0646       theta = atan(1 / cotgth_);
0647     } else {
0648       theta = TMath::Pi() / 2;
0649     }
0650     if (theta < 0)
0651       theta += TMath::Pi();
0652     eta = -log(tan(theta / 2));
0653   }
0654   double cotgth_;
0655   TRandom* gRandom_;
0656 };
0657 inline smearFunctionBase::~smearFunctionBase() {}  // defined even though it's pure virtual; should be faster this way.
0658 
0659 // No smearing
0660 // -----------
0661 class smearFunctionType0 : public smearFunctionBase {
0662 public:
0663   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {}
0664 };
0665 // The 3 parameters of smearType1 are: pt dependence of pt smear, phi smear and
0666 // cotgtheta smear.
0667 class smearFunctionType1 : public smearFunctionBase {
0668 public:
0669   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {
0670     pt = pt * (1.0 + y[0] * parSmear[0] * pt);
0671     phi = phi * (1.0 + y[1] * parSmear[1]);
0672     double tmp = 2 * atan(exp(-eta));
0673     cotgth_ = cos(tmp) / sin(tmp) * (1.0 + y[2] * parSmear[2]);
0674     smearEta(eta);
0675   }
0676 };
0677 
0678 class smearFunctionType2 : public smearFunctionBase {
0679 public:
0680   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {
0681     pt = pt * (1.0 + y[0] * parSmear[0] * pt + y[1] * parSmear[1] * std::fabs(eta));
0682     phi = phi * (1.0 + y[2] * parSmear[2]);
0683     double tmp = 2 * atan(exp(-eta));
0684     cotgth_ = cos(tmp) / sin(tmp) * (1.0 + y[3] * parSmear[3]);
0685     smearEta(eta);
0686   }
0687 };
0688 
0689 class smearFunctionType3 : public smearFunctionBase {
0690 public:
0691   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {
0692     pt = pt * (1.0 + y[0] * parSmear[0] * pt + y[1] * parSmear[1] * std::fabs(eta));
0693     phi = phi * (1.0 + y[2] * parSmear[2]);
0694     double tmp = 2 * atan(exp(-eta));
0695     cotgth_ = cos(tmp) / sin(tmp) * (1.0 + y[3] * parSmear[3] + y[4] * parSmear[4] * std::fabs(eta));
0696     smearEta(eta);
0697   }
0698 };
0699 // The six parameters of SmearType=4 are respectively:
0700 // Pt dep. of Pt res., |eta| dep. of Pt res., Phi res., |eta| res.,
0701 // |eta| dep. of |eta| res., Pt^2 dep. of Pt res.
0702 class smearFunctionType4 : public smearFunctionBase {
0703 public:
0704   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {
0705     pt = pt * (1.0 + y[0] * parSmear[0] * pt + y[1] * parSmear[1] * std::fabs(eta) + y[5] * parSmear[5] * pow(pt, 2));
0706     phi = phi * (1.0 + y[2] * parSmear[2]);
0707     double tmp = 2 * atan(exp(-eta));
0708     cotgth_ = cos(tmp) / sin(tmp) * (1.0 + y[3] * parSmear[3] + y[4] * parSmear[4] * std::fabs(eta));
0709     smearEta(eta);
0710   }
0711 };
0712 
0713 class smearFunctionType5 : public smearFunctionBase {
0714 public:
0715   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {
0716     pt = pt * (1.0 + y[0] * parSmear[0] * pt + y[1] * parSmear[1] * std::fabs(eta) + y[5] * parSmear[5] * pow(pt, 2));
0717     phi = phi * (1.0 + y[2] * parSmear[2] + y[6] * parSmear[6] * pt);
0718     double tmp = 2 * atan(exp(-eta));
0719     cotgth_ = cos(tmp) / sin(tmp) * (1.0 + y[3] * parSmear[3] + y[4] * parSmear[4] * std::fabs(eta));
0720     smearEta(eta);
0721   }
0722 };
0723 
0724 //Smearing for MC correction based on the resolution function Type 15 for misaligned MC
0725 class smearFunctionType6 : public smearFunctionBase {
0726 public:
0727   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {
0728     double sigmaSmear = 0;
0729     double sigmaPtAl = 0;
0730     double sigmaPtMisal = 0;
0731     double ptPart = parSmear[0] + parSmear[1] * 1 / pt + pt * parSmear[2];
0732     double fabsEta = std::fabs(eta);
0733 
0734     sigmaPtAl = parSmear[14] * etaByPoints(eta, parSmear[15]);
0735 
0736     if (std::fabs(eta) <= 1.4) {
0737       sigmaPtMisal = ptPart + parSmear[3] + parSmear[4] * std::fabs(eta) + parSmear[5] * eta * eta;
0738       sigmaSmear = sqrt(std::fabs(pow(sigmaPtMisal, 2) - pow(sigmaPtAl, 2)));
0739       pt = pt * gRandom_->Gaus(1, sigmaSmear);
0740     } else if (eta > 1.4) {  //eta in right endcap
0741       double par =
0742           parSmear[3] + parSmear[4] * 1.4 + parSmear[5] * 1.4 * 1.4 -
0743           (parSmear[6] + parSmear[7] * (1.4 - parSmear[8]) + parSmear[9] * (1.4 - parSmear[8]) * (1.4 - parSmear[8]));
0744       sigmaPtMisal = par + ptPart + parSmear[6] + parSmear[7] * std::fabs((fabsEta - parSmear[8])) +
0745                      parSmear[9] * (fabsEta - parSmear[8]) * (fabsEta - parSmear[8]);
0746       sigmaSmear = sqrt(std::fabs(pow(sigmaPtMisal, 2) - pow(sigmaPtAl, 2)));
0747       pt = pt * gRandom_->Gaus(1, sigmaSmear);
0748     } else {  //eta in left endcap
0749       double par = parSmear[3] + parSmear[4] * 1.4 + parSmear[5] * 1.4 * 1.4 -
0750                    (parSmear[10] + parSmear[11] * (1.4 - parSmear[12]) +
0751                     parSmear[13] * (1.4 - parSmear[12]) * (1.4 - parSmear[12]));
0752       sigmaPtMisal = par + ptPart + parSmear[10] + parSmear[11] * std::fabs((fabsEta - parSmear[12])) +
0753                      parSmear[13] * (fabsEta - parSmear[12]) * (fabsEta - parSmear[12]);
0754       sigmaSmear = sqrt(std::fabs(pow(sigmaPtMisal, 2) - pow(sigmaPtAl, 2)));
0755       pt = pt * gRandom_->Gaus(1, sigmaSmear);
0756     }
0757   }
0758 
0759 protected:
0760   /**
0761    * This is the pt vs eta resolution by points. It uses std::fabs(eta) assuming symmetry.
0762    * The values are derived from 100k events of MuonGun with 5<pt<100 and |eta|<3.
0763    */
0764   double etaByPoints(const double& inEta, const double& border) {
0765     Double_t eta = std::fabs(inEta);
0766     if (0. <= eta && eta <= 0.2)
0767       return 0.00942984;
0768     else if (0.2 < eta && eta <= 0.4)
0769       return 0.0104489;
0770     else if (0.4 < eta && eta <= 0.6)
0771       return 0.0110521;
0772     else if (0.6 < eta && eta <= 0.8)
0773       return 0.0117338;
0774     else if (0.8 < eta && eta <= 1.0)
0775       return 0.0138142;
0776     else if (1.0 < eta && eta <= 1.2)
0777       return 0.0165826;
0778     else if (1.2 < eta && eta <= 1.4)
0779       return 0.0183663;
0780     else if (1.4 < eta && eta <= 1.6)
0781       return 0.0169904;
0782     else if (1.6 < eta && eta <= 1.8)
0783       return 0.0173289;
0784     else if (1.8 < eta && eta <= 2.0)
0785       return 0.0205821;
0786     else if (2.0 < eta && eta <= 2.2)
0787       return 0.0250032;
0788     else if (2.2 < eta && eta <= 2.4)
0789       return 0.0339477;
0790     // ATTENTION: This point has a big error and it is very displaced from the rest of the distribution.
0791     else if (2.4 < eta && eta <= 2.6)
0792       return border;
0793     return (0.);
0794   }
0795 };
0796 
0797 class smearFunctionType7 : public smearFunctionBase {
0798 public:
0799   void smear(double& pt, double& eta, double& phi, const double* y, const std::vector<double>& parSmear) override {
0800     double sigmaSquared = sigmaPtDiff.squaredDiff(eta);
0801     TF1 G("G", "[0]*exp(-0.5*pow(x,2)/[1])", -5., 5.);
0802     double norm = 1 / (sqrt(2 * TMath::Pi() * sigmaSquared));
0803     G.SetParameter(0, norm);
0804     G.SetParameter(1, sigmaSquared);
0805     pt = pt * (1 - G.GetRandom());
0806   }
0807   SigmaPtDiff sigmaPtDiff;
0808 };
0809 
0810 /// Service to build the smearing functor corresponding to the passed identifier
0811 smearFunctionBase* smearFunctionService(const int identifier);
0812 
0813 // // Defined globally...
0814 // static smearFunctionBase * smearFunctionArray[] = {
0815 //   new smearFunctionType0,
0816 //   new smearFunctionType1,
0817 //   new smearFunctionType2,
0818 //   new smearFunctionType3,
0819 //   new smearFunctionType4,
0820 //   new smearFunctionType5
0821 // };
0822 
0823 /**
0824  * Resolution functions. <br>
0825  * Need to use templates to make it work with both array and std::vector<double>.
0826  */
0827 template <class T>
0828 class resolutionFunctionBase {
0829 public:
0830   virtual double sigmaPt(const double& pt, const double& eta, const T& parval) = 0;
0831   virtual double sigmaPtError(const double& pt, const double& eta, const T& parval, const T& parError) { return 0.; }
0832   virtual double sigmaPhi(const double& pt, const double& eta, const T& parval) = 0;
0833   virtual double sigmaCotgTh(const double& pt, const double& eta, const T& parval) = 0;
0834   virtual double covPt1Pt2(
0835       const double& pt1, const double& eta1, const double& pt2, const double& eta2, const T& parval) {
0836     return 0.;
0837   }
0838   resolutionFunctionBase() {}
0839   virtual ~resolutionFunctionBase() = 0;
0840   /// This method is used to differentiate parameters among the different functions
0841   virtual void setParameters(double* Start,
0842                              double* Step,
0843                              double* Mini,
0844                              double* Maxi,
0845                              int* ind,
0846                              TString* parname,
0847                              const T& parResol,
0848                              const std::vector<int>& parResolOrder,
0849                              const int muonType){};
0850   virtual void setParameters(double* Start,
0851                              double* Step,
0852                              double* Mini,
0853                              double* Maxi,
0854                              int* ind,
0855                              TString* parname,
0856                              const T& parResol,
0857                              const std::vector<int>& parResolOrder,
0858                              const std::vector<double>& parStep,
0859                              const std::vector<double>& parMin,
0860                              const std::vector<double>& parMax,
0861                              const int muonType) {
0862     std::cout << "The method setParameters must be implemented for this resolution function" << std::endl;
0863     exit(1);
0864   }
0865   virtual int parNum() const { return parNum_; }
0866 
0867 protected:
0868   int parNum_;
0869   /// This method sets the parameters
0870   virtual void setPar(double* Start,
0871                       double* Step,
0872                       double* Mini,
0873                       double* Maxi,
0874                       int* ind,
0875                       TString* parname,
0876                       const T& parResol,
0877                       const std::vector<int>& parResolOrder,
0878                       double* thisStep,
0879                       double* thisMini,
0880                       double* thisMaxi,
0881                       TString* thisParName) {
0882     for (int iPar = 0; iPar < this->parNum_; ++iPar) {
0883       Start[iPar] = parResol[iPar];
0884       Step[iPar] = thisStep[iPar];
0885       Mini[iPar] = thisMini[iPar];
0886       Maxi[iPar] = thisMaxi[iPar];
0887       ind[iPar] = parResolOrder[iPar];
0888       parname[iPar] = thisParName[iPar];
0889     }
0890   }
0891   virtual void setPar(double* Start,
0892                       double* Step,
0893                       double* Mini,
0894                       double* Maxi,
0895                       int* ind,
0896                       TString* parname,
0897                       const T& parResol,
0898                       const std::vector<int>& parResolOrder,
0899                       const std::vector<ParameterSet>& parSet) {
0900     if (int(parSet.size()) != this->parNum_) {
0901       std::cout << "Error: wrong number of parameter initializations = " << parSet.size()
0902                 << ". Number of parameters is " << this->parNum_ << std::endl;
0903       exit(1);
0904     }
0905     for (int iPar = 0; iPar < this->parNum_; ++iPar) {
0906       Start[iPar] = parResol[iPar];
0907       Step[iPar] = parSet[iPar].step;
0908       Mini[iPar] = parSet[iPar].mini;
0909       Maxi[iPar] = parSet[iPar].maxi;
0910       ind[iPar] = parResolOrder[iPar];
0911       parname[iPar] = parSet[iPar].name;
0912     }
0913   }
0914 };
0915 template <class T>
0916 inline resolutionFunctionBase<T>::~resolutionFunctionBase() {
0917 }  // defined even though it's pure virtual; should be faster this way.
0918 
0919 // SLIMMED VERSION
0920 // The full set of resolutionFunction developed in the past can be found at
0921 // https://github.com/cms-sw/cmssw/blob/CMSSW_5_3_X/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h#L3082-L5052
0922 
0923 // null
0924 // --------
0925 template <class T>
0926 class resolutionFunctionType0 : public resolutionFunctionBase<T> {
0927 public:
0928   resolutionFunctionType0() {
0929     // One of the two is required. This follows from when templates are used by the compiler and the names lookup rules in c++.
0930     this->parNum_ = 0;
0931   }
0932   double sigmaPt(const double& pt, const double& eta, const T& parval) override { return 0; }
0933   double sigmaCotgTh(const double& pt, const double& eta, const T& parval) override { return 0; }
0934   double sigmaPhi(const double& pt, const double& eta, const T& parval) override { return 0.; }
0935   // derivatives ---------------
0936   double sigmaPtError(const double& pt, const double& eta, const T& parval, const T& parError) override { return 0; }
0937   void setParameters(double* Start,
0938                      double* Step,
0939                      double* Mini,
0940                      double* Maxi,
0941                      int* ind,
0942                      TString* parname,
0943                      const T& parResol,
0944                      const std::vector<int>& parResolOrder,
0945                      const int muonType) override {}
0946   void setParameters(double* Start,
0947                      double* Step,
0948                      double* Mini,
0949                      double* Maxi,
0950                      int* ind,
0951                      TString* parname,
0952                      const T& parResol,
0953                      const std::vector<int>& parResolOrder,
0954                      const std::vector<double>& parStep,
0955                      const std::vector<double>& parMin,
0956                      const std::vector<double>& parMax,
0957                      const int muonType) override {}
0958 };
0959 
0960 // Binned in eta to fit the Z (parametrization as linear sum)
0961 template <class T>
0962 class resolutionFunctionType45 : public resolutionFunctionBase<T> {
0963 public:
0964   int etaBin(const double& eta) {
0965     // 12 bins from -2.4 to 2.4
0966     // std::cout << "for eta = " << eta << ", bin = " << bin << std::endl;
0967     if (eta < -2.0)
0968       return 1;
0969     if (eta < -1.8)
0970       return 2;
0971     if (eta < -1.6)
0972       return 3;
0973     if (eta < -1.2)
0974       return 4;
0975     if (eta < -0.8)
0976       return 5;
0977     if (eta < 0.)
0978       return 6;
0979     if (eta < 0.8)
0980       return 7;
0981     if (eta < 1.2)
0982       return 8;
0983     if (eta < 1.6)
0984       return 9;
0985     if (eta < 1.8)
0986       return 10;
0987     if (eta < 2.0)
0988       return 11;
0989     return 12;
0990   }
0991 
0992   resolutionFunctionType45() { this->parNum_ = 13; }
0993 
0994   double sigmaPt(const double& pt, const double& eta, const T& parval) override {
0995     return (parval[0] * pt + parval[etaBin(eta)]);
0996   }
0997   double sigmaCotgTh(const double& pt, const double& eta, const T& parval) override { return 0; }
0998   double sigmaPhi(const double& pt, const double& eta, const T& parval) override { return 0.; }
0999 
1000   // derivatives ---------------
1001   double sigmaPtError(const double& pt, const double& eta, const T& parval, const T& parError) override {
1002     // Use the etaBin function to select the right bin for the parameter
1003     return sqrt(pow(pt * parError[0], 2) + pow(parError[etaBin(eta)], 2));
1004   }
1005 
1006   void setParameters(double* Start,
1007                      double* Step,
1008                      double* Mini,
1009                      double* Maxi,
1010                      int* ind,
1011                      TString* parname,
1012                      const T& parResol,
1013                      const std::vector<int>& parResolOrder,
1014                      const int muonType) override {
1015     std::vector<ParameterSet> parSet(this->parNum_);
1016     // name, step, mini, maxi
1017     parSet[0] = ParameterSet("Pt res. sc.", 0.002, -0.1, 0.1);
1018     parSet[1] = ParameterSet("eta bin 1", 0.00002, -0.01, 0.01);
1019     parSet[2] = ParameterSet("eta bin 2", 0.00002, -0.01, 0.01);
1020     parSet[3] = ParameterSet("eta bin 3", 0.00002, -0.01, 0.01);
1021     parSet[4] = ParameterSet("eta bin 4", 0.00002, -0.01, 0.01);
1022     parSet[5] = ParameterSet("eta bin 5", 0.00002, -0.01, 0.01);
1023     parSet[6] = ParameterSet("eta bin 6", 0.00002, -0.01, 0.01);
1024     parSet[7] = ParameterSet("eta bin 7", 0.00002, -0.01, 0.01);
1025     parSet[8] = ParameterSet("eta bin 8", 0.00002, -0.01, 0.01);
1026     parSet[9] = ParameterSet("eta bin 9", 0.00002, -0.01, 0.01);
1027     parSet[10] = ParameterSet("eta bin 10", 0.00002, -0.01, 0.01);
1028     parSet[11] = ParameterSet("eta bin 11", 0.00002, -0.01, 0.01);
1029     parSet[12] = ParameterSet("eta bin 12", 0.00002, -0.01, 0.01);
1030 
1031     std::cout << "setting parameters" << std::endl;
1032     this->setPar(Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parSet);
1033   }
1034 
1035   void setParameters(double* Start,
1036                      double* Step,
1037                      double* Mini,
1038                      double* Maxi,
1039                      int* ind,
1040                      TString* parname,
1041                      const T& parResol,
1042                      const std::vector<int>& parResolOrder,
1043                      const std::vector<double>& parStep,
1044                      const std::vector<double>& parMin,
1045                      const std::vector<double>& parMax,
1046                      const int muonType) override {
1047     if ((int(parStep.size()) != this->parNum_) || (int(parMin.size()) != this->parNum_) ||
1048         (int(parMax.size()) != this->parNum_)) {
1049       std::cout << "Error: par step or min or max do not match with number of parameters" << std::endl;
1050       std::cout << "parNum = " << this->parNum_ << std::endl;
1051       std::cout << "parStep.size() = " << parStep.size() << std::endl;
1052       std::cout << "parMin.size() = " << parMin.size() << std::endl;
1053       std::cout << "parMax.size() = " << parMax.size() << std::endl;
1054       exit(1);
1055     }
1056     std::vector<ParameterSet> parSet(this->parNum_);
1057     // name, step, mini, maxi
1058     parSet[0] = ParameterSet("Pt res. sc.", parStep[0], parMin[0], parMax[0]);
1059     parSet[1] = ParameterSet("eta bin 1", parStep[1], parMin[1], parMax[1]);
1060     parSet[2] = ParameterSet("eta bin 2", parStep[2], parMin[2], parMax[2]);
1061     parSet[3] = ParameterSet("eta bin 3", parStep[3], parMin[3], parMax[3]);
1062     parSet[4] = ParameterSet("eta bin 4", parStep[4], parMin[4], parMax[4]);
1063     parSet[5] = ParameterSet("eta bin 5", parStep[5], parMin[5], parMax[5]);
1064     parSet[6] = ParameterSet("eta bin 6", parStep[6], parMin[6], parMax[6]);
1065     parSet[7] = ParameterSet("eta bin 7", parStep[7], parMin[7], parMax[7]);
1066     parSet[8] = ParameterSet("eta bin 8", parStep[8], parMin[8], parMax[8]);
1067     parSet[9] = ParameterSet("eta bin 9", parStep[9], parMin[9], parMax[9]);
1068     parSet[10] = ParameterSet("eta bin 10", parStep[10], parMin[10], parMax[10]);
1069     parSet[11] = ParameterSet("eta bin 11", parStep[11], parMin[11], parMax[11]);
1070     parSet[12] = ParameterSet("eta bin 12", parStep[12], parMin[12], parMax[12]);
1071 
1072     std::cout << "setting parameters" << std::endl;
1073     for (int i = 0; i < this->parNum_; ++i) {
1074       std::cout << "parStep[" << i << "] = " << parStep[i] << ", parMin[" << i << "] = " << parMin[i] << ", parMax["
1075                 << i << "] = " << parMin[i] << std::endl;
1076     }
1077     this->setPar(Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parSet);
1078   }
1079 };
1080 
1081 // Binned in eta to fit the Z (parametrization as sum in quadrature)
1082 template <class T>
1083 class resolutionFunctionType46 : public resolutionFunctionBase<T> {
1084 public:
1085   int etaBin(const double& eta) {
1086     // 12 bins from -2.4 to 2.4
1087     // std::cout << "for eta = " << eta << ", bin = " << bin << std::endl;
1088     if (eta < -2.0)
1089       return 1;
1090     if (eta < -1.8)
1091       return 2;
1092     if (eta < -1.6)
1093       return 3;
1094     if (eta < -1.2)
1095       return 4;
1096     if (eta < -0.8)
1097       return 5;
1098     if (eta < 0.)
1099       return 6;
1100     if (eta < 0.8)
1101       return 7;
1102     if (eta < 1.2)
1103       return 8;
1104     if (eta < 1.6)
1105       return 9;
1106     if (eta < 1.8)
1107       return 10;
1108     if (eta < 2.0)
1109       return 11;
1110     return 12;
1111   }
1112 
1113   resolutionFunctionType46() { this->parNum_ = 13; }
1114 
1115   double sigmaPt(const double& pt, const double& eta, const T& parval) override {
1116     return sqrt(pow(parval[0] * pt, 2) + pow(parval[etaBin(eta)], 2));
1117   }
1118   double sigmaCotgTh(const double& pt, const double& eta, const T& parval) override { return 0; }
1119   double sigmaPhi(const double& pt, const double& eta, const T& parval) override { return 0.; }
1120 
1121   // derivatives ---------------
1122   double sigmaPtError(const double& pt, const double& eta, const T& parval, const T& parError) override {
1123     // Use the etaByBin function to select the right bin for the parameter
1124     double r = sqrt(pow(parval[0] * pt, 2) + pow(parval[etaBin(eta)], 2));
1125     return sqrt(pow(pt * pt * parval[0] * parError[0], 2) + pow(parval[etaBin(eta)] * parError[etaBin(eta)], 2)) / r;
1126   }
1127 
1128   void setParameters(double* Start,
1129                      double* Step,
1130                      double* Mini,
1131                      double* Maxi,
1132                      int* ind,
1133                      TString* parname,
1134                      const T& parResol,
1135                      const std::vector<int>& parResolOrder,
1136                      const int muonType) override {
1137     std::vector<ParameterSet> parSet(this->parNum_);
1138     // name, step, mini, maxi
1139     parSet[0] = ParameterSet("Pt res. sc.", 0.0002, 0., 0.1);
1140     parSet[1] = ParameterSet("eta bin 1", 0.00002, 0., 0.01);
1141     parSet[2] = ParameterSet("eta bin 2", 0.00002, 0., 0.01);
1142     parSet[3] = ParameterSet("eta bin 3", 0.00002, 0., 0.01);
1143     parSet[4] = ParameterSet("eta bin 4", 0.00002, 0., 0.01);
1144     parSet[5] = ParameterSet("eta bin 5", 0.00002, 0., 0.01);
1145     parSet[6] = ParameterSet("eta bin 6", 0.00002, 0., 0.01);
1146     parSet[7] = ParameterSet("eta bin 7", 0.00002, 0., 0.01);
1147     parSet[8] = ParameterSet("eta bin 8", 0.00002, 0., 0.01);
1148     parSet[9] = ParameterSet("eta bin 9", 0.00002, 0., 0.01);
1149     parSet[10] = ParameterSet("eta bin 10", 0.00002, 0., 0.01);
1150     parSet[11] = ParameterSet("eta bin 11", 0.00002, 0., 0.01);
1151     parSet[12] = ParameterSet("eta bin 12", 0.00002, 0., 0.01);
1152 
1153     std::cout << "setting parameters" << std::endl;
1154     this->setPar(Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parSet);
1155   }
1156 
1157   void setParameters(double* Start,
1158                      double* Step,
1159                      double* Mini,
1160                      double* Maxi,
1161                      int* ind,
1162                      TString* parname,
1163                      const T& parResol,
1164                      const std::vector<int>& parResolOrder,
1165                      const std::vector<double>& parStep,
1166                      const std::vector<double>& parMin,
1167                      const std::vector<double>& parMax,
1168                      const int muonType) override {
1169     if ((int(parStep.size()) != this->parNum_) || (int(parMin.size()) != this->parNum_) ||
1170         (int(parMax.size()) != this->parNum_)) {
1171       std::cout << "Error: par step or min or max do not match with number of parameters" << std::endl;
1172       std::cout << "parNum = " << this->parNum_ << std::endl;
1173       std::cout << "parStep.size() = " << parStep.size() << std::endl;
1174       std::cout << "parMin.size() = " << parMin.size() << std::endl;
1175       std::cout << "parMax.size() = " << parMax.size() << std::endl;
1176       exit(1);
1177     }
1178     std::vector<ParameterSet> parSet(this->parNum_);
1179     // name, step, mini, maxi
1180     parSet[0] = ParameterSet("Pt res. sc.", parStep[0], parMin[0], parMax[0]);
1181     parSet[1] = ParameterSet("eta bin 1", parStep[1], parMin[1], parMax[1]);
1182     parSet[2] = ParameterSet("eta bin 2", parStep[2], parMin[2], parMax[2]);
1183     parSet[3] = ParameterSet("eta bin 3", parStep[3], parMin[3], parMax[3]);
1184     parSet[4] = ParameterSet("eta bin 4", parStep[4], parMin[4], parMax[4]);
1185     parSet[5] = ParameterSet("eta bin 5", parStep[5], parMin[5], parMax[5]);
1186     parSet[6] = ParameterSet("eta bin 6", parStep[6], parMin[6], parMax[6]);
1187     parSet[7] = ParameterSet("eta bin 7", parStep[7], parMin[7], parMax[7]);
1188     parSet[8] = ParameterSet("eta bin 8", parStep[8], parMin[8], parMax[8]);
1189     parSet[9] = ParameterSet("eta bin 9", parStep[9], parMin[9], parMax[9]);
1190     parSet[10] = ParameterSet("eta bin 10", parStep[10], parMin[10], parMax[10]);
1191     parSet[11] = ParameterSet("eta bin 11", parStep[11], parMin[11], parMax[11]);
1192     parSet[12] = ParameterSet("eta bin 12", parStep[12], parMin[12], parMax[12]);
1193 
1194     std::cout << "setting parameters" << std::endl;
1195     for (int i = 0; i < this->parNum_; ++i) {
1196       std::cout << "parStep[" << i << "] = " << parStep[i] << ", parMin[" << i << "] = " << parMin[i] << ", parMax["
1197                 << i << "] = " << parMin[i] << std::endl;
1198     }
1199     this->setPar(Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parSet);
1200   }
1201 };
1202 
1203 // Binned in eta to fit the Z (parametrization as sum in quadrature) and including an overall covariance
1204 template <class T>
1205 class resolutionFunctionType47 : public resolutionFunctionBase<T> {
1206 public:
1207   int etaBin(const double& eta) {
1208     // 12 bins from -2.4 to 2.4
1209     // std::cout << "for eta = " << eta << ", bin = " << bin << std::endl;
1210     if (eta < -2.0)
1211       return 1;
1212     if (eta < -1.8)
1213       return 2;
1214     if (eta < -1.6)
1215       return 3;
1216     if (eta < -1.2)
1217       return 4;
1218     if (eta < -0.8)
1219       return 5;
1220     if (eta < 0.)
1221       return 6;
1222     if (eta < 0.8)
1223       return 7;
1224     if (eta < 1.2)
1225       return 8;
1226     if (eta < 1.6)
1227       return 9;
1228     if (eta < 1.8)
1229       return 10;
1230     if (eta < 2.0)
1231       return 11;
1232     return 12;
1233   }
1234 
1235   resolutionFunctionType47() { this->parNum_ = 14; }
1236 
1237   double sigmaPt(const double& pt, const double& eta, const T& parval) override {
1238     return sqrt(pow(parval[0] * pt, 2) + pow(parval[etaBin(eta)], 2) + pow(parval[13], 2));
1239   }
1240   double sigmaCotgTh(const double& pt, const double& eta, const T& parval) override { return 0; }
1241   double sigmaPhi(const double& pt, const double& eta, const T& parval) override { return 0.; }
1242 
1243   double covPt1Pt2(
1244       const double& pt1, const double& eta1, const double& pt2, const double& eta2, const T& parval) override {
1245     return parval[14];
1246   }
1247 
1248   // derivatives ---------------
1249   double sigmaPtError(const double& pt, const double& eta, const T& parval, const T& parError) override {
1250     // Use the etaByBin function to select the right bin for the parameter
1251     double r = sqrt(pow(parval[0] * pt, 2) + pow(parval[etaBin(eta)], 2) + pow(parval[13], 2));
1252     return sqrt(pow(pt * pt * parval[0] * parError[0], 2) + pow(parval[etaBin(eta)] * parError[etaBin(eta)], 2) +
1253                 pow(parval[13] * parError[13], 2)) /
1254            r;
1255   }
1256 
1257   void setParameters(double* Start,
1258                      double* Step,
1259                      double* Mini,
1260                      double* Maxi,
1261                      int* ind,
1262                      TString* parname,
1263                      const T& parResol,
1264                      const std::vector<int>& parResolOrder,
1265                      const int muonType) override {
1266     std::vector<ParameterSet> parSet(this->parNum_);
1267     // name, step, mini, maxi
1268     parSet[0] = ParameterSet("Pt res. sc.", 0.0002, 0., 0.1);
1269     parSet[1] = ParameterSet("eta bin 1", 0.00002, 0., 0.01);
1270     parSet[2] = ParameterSet("eta bin 2", 0.00002, 0., 0.01);
1271     parSet[3] = ParameterSet("eta bin 3", 0.00002, 0., 0.01);
1272     parSet[4] = ParameterSet("eta bin 4", 0.00002, 0., 0.01);
1273     parSet[5] = ParameterSet("eta bin 5", 0.00002, 0., 0.01);
1274     parSet[6] = ParameterSet("eta bin 6", 0.00002, 0., 0.01);
1275     parSet[7] = ParameterSet("eta bin 7", 0.00002, 0., 0.01);
1276     parSet[8] = ParameterSet("eta bin 8", 0.00002, 0., 0.01);
1277     parSet[9] = ParameterSet("eta bin 9", 0.00002, 0., 0.01);
1278     parSet[10] = ParameterSet("eta bin 10", 0.00002, 0., 0.01);
1279     parSet[11] = ParameterSet("eta bin 11", 0.00002, 0., 0.01);
1280     parSet[12] = ParameterSet("eta bin 12", 0.00002, 0., 0.01);
1281     parSet[13] = ParameterSet("cov(pt1,pt2)", 0.00002, 0., 0.01);
1282 
1283     std::cout << "setting parameters" << std::endl;
1284     this->setPar(Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parSet);
1285   }
1286 
1287   void setParameters(double* Start,
1288                      double* Step,
1289                      double* Mini,
1290                      double* Maxi,
1291                      int* ind,
1292                      TString* parname,
1293                      const T& parResol,
1294                      const std::vector<int>& parResolOrder,
1295                      const std::vector<double>& parStep,
1296                      const std::vector<double>& parMin,
1297                      const std::vector<double>& parMax,
1298                      const int muonType) override {
1299     if ((int(parStep.size()) != this->parNum_) || (int(parMin.size()) != this->parNum_) ||
1300         (int(parMax.size()) != this->parNum_)) {
1301       std::cout << "Error: par step or min or max do not match with number of parameters" << std::endl;
1302       std::cout << "parNum = " << this->parNum_ << std::endl;
1303       std::cout << "parStep.size() = " << parStep.size() << std::endl;
1304       std::cout << "parMin.size() = " << parMin.size() << std::endl;
1305       std::cout << "parMax.size() = " << parMax.size() << std::endl;
1306       exit(1);
1307     }
1308     std::vector<ParameterSet> parSet(this->parNum_);
1309     // name, step, mini, maxi
1310     parSet[0] = ParameterSet("Pt res. sc.", parStep[0], parMin[0], parMax[0]);
1311     parSet[1] = ParameterSet("eta bin 1", parStep[1], parMin[1], parMax[1]);
1312     parSet[2] = ParameterSet("eta bin 2", parStep[2], parMin[2], parMax[2]);
1313     parSet[3] = ParameterSet("eta bin 3", parStep[3], parMin[3], parMax[3]);
1314     parSet[4] = ParameterSet("eta bin 4", parStep[4], parMin[4], parMax[4]);
1315     parSet[5] = ParameterSet("eta bin 5", parStep[5], parMin[5], parMax[5]);
1316     parSet[6] = ParameterSet("eta bin 6", parStep[6], parMin[6], parMax[6]);
1317     parSet[7] = ParameterSet("eta bin 7", parStep[7], parMin[7], parMax[7]);
1318     parSet[8] = ParameterSet("eta bin 8", parStep[8], parMin[8], parMax[8]);
1319     parSet[9] = ParameterSet("eta bin 9", parStep[9], parMin[9], parMax[9]);
1320     parSet[10] = ParameterSet("eta bin 10", parStep[10], parMin[10], parMax[10]);
1321     parSet[11] = ParameterSet("eta bin 11", parStep[11], parMin[11], parMax[11]);
1322     parSet[12] = ParameterSet("eta bin 12", parStep[12], parMin[12], parMax[12]);
1323     parSet[13] = ParameterSet("cov(pt1,pt2)", parStep[13], parMin[13], parMax[13]);
1324 
1325     std::cout << "setting parameters" << std::endl;
1326     for (int i = 0; i < this->parNum_; ++i) {
1327       std::cout << "parStep[" << i << "] = " << parStep[i] << ", parMin[" << i << "] = " << parMin[i] << ", parMax["
1328                 << i << "] = " << parMin[i] << std::endl;
1329     }
1330     this->setPar(Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parSet);
1331   }
1332 };
1333 
1334 // ------------ ATTENTION ----------- //
1335 // Other functions are not in for now //
1336 // ---------------------------------- //
1337 
1338 /// Service to build the resolution functor corresponding to the passed identifier
1339 resolutionFunctionBase<double*>* resolutionFunctionService(const int identifier);
1340 
1341 /// Service to build the resolution functor corresponding to the passed identifier when receiving a std::vector<double>
1342 resolutionFunctionBase<std::vector<double> >* resolutionFunctionVecService(const int identifier);
1343 
1344 /**
1345  * Background functors. <br>
1346  * MuScleFit uses different background functions for each resonance. This is done because the
1347  * background shape can change from a region to another so that it is not well described just one of the following functions.
1348  * Since we are only interested in getting the correct shape and fraction in the resonance region we can split it in several
1349  * parts and fit them separately. <br>
1350  * <br>
1351  * When fitting the background function: <br>
1352  * - we define three regions: <br>
1353  * -- Psis region: centered in the mean value of J/Psi and Psi2S masses. <br>
1354  * -- Upsilons region: centered in the mean value of all the Upsilon masses. <br>
1355  * -- Z region. <br>
1356  * - In this case we thus have only three background functions, that is three sets of parameters. <br>
1357  * When not fitting the background function: <br>
1358  * - the windows considered are much narrower and centered around each resonance. <br>
1359  * - In this case we compute a rescaled background fraction parameter for each resonance and we use the rest of the parameters
1360  * for the corresponding region. <br>
1361  * ATTENTION: we are assuming that J/Psi and Psi2S have the same bin width. The same assumption is done for the Upsilons. <br>
1362  * This is the case in the present probability root file. If this changes, the background function for the different regions
1363  * must keep it into account. <br>
1364  * <br>
1365  * All this is handled in MuScleFit by a single functor: BackgroundHandler defined in another file in this dir. <br>
1366  * ATTENTION: The BackgroundHandler assumes that the background fraction is always the first parameter (parBgr[0]).
1367  */
1368 class backgroundFunctionBase {
1369 public:
1370   backgroundFunctionBase(const double& lowerLimit, const double& upperLimit)
1371       : lowerLimit_(lowerLimit), upperLimit_(upperLimit) {}
1372   virtual ~backgroundFunctionBase() { delete functionForIntegral_; };
1373   virtual double operator()(const double* parval, const double& mass, const double& eta) const = 0;
1374   virtual double operator()(const double* parval, const double& mass, const double& eta1, const double& eta2) const {
1375     return operator()(parval, mass, eta1);
1376   }
1377   virtual int parNum() const { return parNum_; }
1378   /// This method is used to differentiate parameters among the different functions
1379   virtual void setParameters(double* Start,
1380                              double* Step,
1381                              double* Mini,
1382                              double* Maxi,
1383                              int* ind,
1384                              TString* parname,
1385                              const std::vector<double>::const_iterator& parBgrIt,
1386                              const std::vector<int>::const_iterator& parBgrOrderIt,
1387                              const int muonType) = 0;
1388   virtual TF1* functionForIntegral(const std::vector<double>::const_iterator& parBgrIt) const {
1389     functionForIntegral_ = new FunctionForIntegral(this, parBgrIt);
1390     TF1* backgroundFunctionForIntegral =
1391         new TF1("backgroundFunctionForIntegral", functionForIntegral_, lowerLimit_, upperLimit_, this->parNum_);
1392     return (backgroundFunctionForIntegral);
1393   }
1394   virtual double fracVsEta(const double* parval, const double& eta1, const double& eta2) const { return 1.; }
1395 
1396 protected:
1397   int parNum_;
1398   double lowerLimit_;
1399   double upperLimit_;
1400   /// This method sets the parameters
1401   virtual void setPar(double* Start,
1402                       double* Step,
1403                       double* Mini,
1404                       double* Maxi,
1405                       int* ind,
1406                       TString* parname,
1407                       const std::vector<double>::const_iterator& parBgrIt,
1408                       const std::vector<int>::const_iterator& parBgrOrderIt,
1409                       double* thisStep,
1410                       double* thisMini,
1411                       double* thisMaxi,
1412                       TString* thisParName) {
1413     for (int iPar = 0; iPar < this->parNum_; ++iPar) {
1414       Start[iPar] = *(parBgrIt + iPar);
1415       Step[iPar] = thisStep[iPar];
1416       Mini[iPar] = thisMini[iPar];
1417       Maxi[iPar] = thisMaxi[iPar];
1418       ind[iPar] = *(parBgrOrderIt + iPar);
1419       // EM 2012.05.22 this line is crashing cmsRun (need to be fixed)      parname[iPar] = thisParName[iPar];
1420     }
1421   }
1422   class FunctionForIntegral {
1423   public:
1424     FunctionForIntegral(const backgroundFunctionBase* function, const std::vector<double>::const_iterator& parBgrIt)
1425         : function_(function) {
1426       parval_ = new double[function_->parNum()];
1427       for (int i = 0; i < function_->parNum(); ++i) {
1428         parval_[i] = *(parBgrIt + i);
1429       }
1430     }
1431     ~FunctionForIntegral() { delete parval_; }
1432     double operator()(const double* mass, const double*) const {
1433       // FIXME: this is a gross approximation. The function should be integrated in eta over the sample.
1434       return ((*function_)(parval_, *mass, 0.));
1435     }
1436 
1437   protected:
1438     const backgroundFunctionBase* function_;
1439     double* parval_;
1440   };
1441   mutable FunctionForIntegral* functionForIntegral_;
1442 };
1443 
1444 /// Linear
1445 // -------
1446 class backgroundFunctionType1 : public backgroundFunctionBase {
1447 public:
1448   /**
1449    * Returns the value of the linear function f(M) = 1 + b*M for M < -1/b, 0 otherwise. <br>
1450    * b is chosen to be negative (background decreasing when M increases). <br>
1451    * Note that this form describes only cases with a != 0 (keep in mind that the relative height
1452    * with respect to the signal is controlled by the fraction parameter).
1453    */
1454   backgroundFunctionType1(const double& lowerLimit, const double& upperLimit)
1455       : backgroundFunctionBase(lowerLimit, upperLimit) {
1456     this->parNum_ = 2;
1457   }
1458   double operator()(const double* parval, const double& mass, const double& eta) const override {
1459     double a = 1.;
1460     double b = parval[1];
1461 
1462     double norm = -(a * lowerLimit_ + b * lowerLimit_ * lowerLimit_ / 2.);
1463 
1464     if (-a / b > upperLimit_)
1465       norm += a * upperLimit_ + b * upperLimit_ * upperLimit_ / 2.;
1466     else
1467       norm += -a * a / (2 * b);
1468 
1469     if (mass < -a / b && norm != 0)
1470       return (a + b * mass) / norm;
1471     else
1472       return 0;
1473   }
1474   void setParameters(double* Start,
1475                      double* Step,
1476                      double* Mini,
1477                      double* Maxi,
1478                      int* ind,
1479                      TString* parname,
1480                      const std::vector<double>::const_iterator& parBgrIt,
1481                      const std::vector<int>::const_iterator& parBgrOrderIt,
1482                      const int muonType) override {
1483     double thisStep[] = {0.01, 0.01};
1484     TString thisParName[] = {"Constant", "Linear"};
1485     if (muonType == 1) {
1486       double thisMini[] = {0.0, -300.};
1487       double thisMaxi[] = {1.0, 0.};
1488       this->setPar(
1489           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1490     } else {
1491       double thisMini[] = {0.0, -300.};
1492       double thisMaxi[] = {1.0, 0.};
1493       this->setPar(
1494           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1495     }
1496   }
1497 };
1498 /// Exponential
1499 // ------------
1500 class backgroundFunctionType2 : public backgroundFunctionBase {
1501 public:
1502   /**
1503    * In case of an exponential, we normalize it such that it has integral in any window
1504    * equal to unity, and then, when adding together all the resonances, one gets a meaningful
1505    * result for the overall background fraction.
1506    */
1507   backgroundFunctionType2(const double& lowerLimit, const double& upperLimit)
1508       : backgroundFunctionBase(lowerLimit, upperLimit) {
1509     this->parNum_ = 2;
1510   }
1511   double operator()(const double* parval, const double& mass, const double& eta) const override {
1512     double Bgrp2 = parval[1];
1513     double norm = -(exp(-Bgrp2 * upperLimit_) - exp(-Bgrp2 * lowerLimit_)) / Bgrp2;
1514     if (norm != 0)
1515       return exp(-Bgrp2 * mass) / norm;
1516     else
1517       return 0.;
1518   }
1519   void setParameters(double* Start,
1520                      double* Step,
1521                      double* Mini,
1522                      double* Maxi,
1523                      int* ind,
1524                      TString* parname,
1525                      const std::vector<double>::const_iterator& parBgrIt,
1526                      const std::vector<int>::const_iterator& parBgrOrderIt,
1527                      const int muonType) override {
1528     double thisStep[] = {0.01, 0.01};
1529     TString thisParName[] = {"Bgr fraction", "Bgr slope"};
1530     if (muonType == 1) {
1531       double thisMini[] = {0.0, 0.};
1532       double thisMaxi[] = {1.0, 10.};
1533       this->setPar(
1534           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1535     } else {
1536       double thisMini[] = {0.0, 0.};
1537       double thisMaxi[] = {1.0, 10.};
1538       this->setPar(
1539           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1540     }
1541   }
1542 
1543   // virtual double fracVsEta(const double * parval, const double & resEta) const
1544   // {
1545   //   // return( 0.6120 - 0.0225*eta*eta );
1546   //   return( 1. - 0.0225*resEta*resEta ); // so that a = 1 for eta = 0.
1547   // }
1548 };
1549 
1550 ///// Constant + Exponential
1551 //// -------------------------------------------------------------------------------- //
1552 //// ATTENTION: TODO: the normalization must be adapted to the asymmetric mass window //
1553 //// -------------------------------------------------------------------------------- //
1554 //
1555 //class backgroundFunctionType3 : public backgroundFunctionBase {
1556 // public:
1557 //  // pass parval[shift]
1558 //  backgroundFunctionType3(const double & lowerLimit, const double & upperLimit) :
1559 //    backgroundFunctionBase(lowerLimit, upperLimit)
1560 //    { this->parNum_ = 3; }
1561 //  virtual double operator()( const double * parval, const int resTotNum, const int ires, const bool * resConsidered,
1562 //                             const double * ResMass, const double ResHalfWidth[], const int MuonType, const double & mass, const int nbins ) {
1563 //    double PB = 0.;
1564 //    double Bgrp2 = parval[1];
1565 //    double Bgrp3 = parval[2];
1566 //    for (int ires=0; ires<resTotNum; ires++) {
1567 //      // In this case, by integrating between A and B, we get for f=exp(a-bx)+k:
1568 //      // INT = exp(a)/b*(exp(-bA)-exp(-bB))+k*(B-A) so our function, which in 1000 bins between A and B
1569 //      // gets a total of 1, is f = (exp(a-bx)+k)*(B-A)/nbins / (INT)
1570 //      // ----------------------------------------------------------------------------------------------
1571 //      if (resConsidered[ires]) {
1572 //  if (exp(-Bgrp2*(ResMass[ires]-ResHalfWidth[ires]))-exp(-Bgrp2*(ResMass[ires]+ResHalfWidth[ires]))>0) {
1573 //    PB += (exp(-Bgrp2*mass)+Bgrp3) *
1574 //      2*ResHalfWidth[ires]/(double)nbins /
1575 //      ( (exp(-Bgrp2*(ResMass[ires]-ResHalfWidth[ires]))-exp(-Bgrp2*(ResMass[ires]+ResHalfWidth[ires])))/
1576 //        Bgrp2 + Bgrp3*2*ResHalfWidth[ires] );
1577 //  } else {
1578 //    std::cout << "Impossible to compute Background probability! - some fix needed - Bgrp2=" << Bgrp2 << std::endl;
1579 //  }
1580 //      }
1581 //    }
1582 //    return PB;
1583 //  }
1584 //  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) {
1585 //    double thisStep[] = {0.1, 0.001, 0.1};
1586 //    TString thisParName[] = {"Bgr fraction", "Bgr slope", "Bgr constant"};
1587 //    if( muonType == 1 ) {
1588 //      double thisMini[] = {0.0, 0.000000001, 0.0};
1589 //      double thisMaxi[] = {1.0, 0.2, 1000};
1590 //      this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
1591 //    } else {
1592 //      double thisMini[] = {0.0, 0.000000001, 0.0};
1593 //      double thisMaxi[] = {1.0, 0.2, 1000};
1594 //      this->setPar( Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName );
1595 //    }
1596 //  }
1597 //  virtual TF1* functionForIntegral(const std::vector<double>::const_iterator & parBgrIt) const {return 0;};
1598 //};
1599 
1600 /// Exponential with eta dependence
1601 // --------------------------------
1602 class backgroundFunctionType4 : public backgroundFunctionBase {
1603 public:
1604   /**
1605    * In case of an exponential, we normalize it such that it has integral in any window
1606    * equal to unity, and then, when adding together all the resonances, one gets a meaningful
1607    * result for the overall background fraction.
1608    */
1609   backgroundFunctionType4(const double& lowerLimit, const double& upperLimit)
1610       : backgroundFunctionBase(lowerLimit, upperLimit) {
1611     this->parNum_ = 4;
1612   }
1613   double operator()(const double* parval, const double& mass, const double& eta) const override {
1614     double Bgrp2 = parval[1] + parval[2] * eta * eta;
1615     double norm = -(exp(-Bgrp2 * upperLimit_) - exp(-Bgrp2 * lowerLimit_)) / Bgrp2;
1616     if (norm != 0)
1617       return exp(-Bgrp2 * mass) / norm;
1618     else
1619       return 0.;
1620   }
1621   void setParameters(double* Start,
1622                      double* Step,
1623                      double* Mini,
1624                      double* Maxi,
1625                      int* ind,
1626                      TString* parname,
1627                      const std::vector<double>::const_iterator& parBgrIt,
1628                      const std::vector<int>::const_iterator& parBgrOrderIt,
1629                      const int muonType) override {
1630     double thisStep[] = {0.01, 0.01, 0.01, 0.01};
1631     TString thisParName[] = {
1632         "Bgr fraction", "Bgr slope", "Bgr slope eta^2 dependence", "background fraction eta dependence"};
1633     if (muonType == 1) {
1634       double thisMini[] = {0.0, 0., 0., -1.};
1635       double thisMaxi[] = {1.0, 10., 10., 1.};
1636       this->setPar(
1637           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1638     } else {
1639       double thisMini[] = {0.0, 0., -1., -1.};
1640       double thisMaxi[] = {1.0, 10., 1., 1.};
1641       this->setPar(
1642           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1643     }
1644   }
1645   /* virtual double fracVsEta(const double * parval, const double & resEta) const */
1646   /* { */
1647   /*   return( 1. - parval[3]*resEta*resEta ); // so that a = 1 for eta = 0. */
1648   /* } */
1649 };
1650 
1651 /// Linear with eta dependence
1652 // ---------------------------
1653 class backgroundFunctionType5 : public backgroundFunctionBase {
1654 public:
1655   /**
1656    * Returns the value of the linear function f(M) = a + b*M for M < -a/b, 0 otherwise. <br>
1657    * Where a = 1 + c*eta*eta and b is chosen to be negative (background decreasing when M increases).
1658    */
1659   backgroundFunctionType5(const double& lowerLimit, const double& upperLimit)
1660       : backgroundFunctionBase(lowerLimit, upperLimit) {
1661     this->parNum_ = 3;
1662   }
1663   double operator()(const double* parval, const double& mass, const double& eta) const override {
1664     double b = parval[1];
1665     // double c = parval[2];
1666     double a = 1 + parval[2] * eta * eta;
1667 
1668     double norm = -(a * lowerLimit_ + b * lowerLimit_ * lowerLimit_ / 2.);
1669 
1670     if (-a / b > upperLimit_)
1671       norm += a * upperLimit_ + b * upperLimit_ * upperLimit_ / 2.;
1672     else
1673       norm += -a * a / (2 * b);
1674 
1675     if (mass < -a / b && norm != 0)
1676       return (a + b * mass) / norm;
1677     else
1678       return 0;
1679   }
1680   void setParameters(double* Start,
1681                      double* Step,
1682                      double* Mini,
1683                      double* Maxi,
1684                      int* ind,
1685                      TString* parname,
1686                      const std::vector<double>::const_iterator& parBgrIt,
1687                      const std::vector<int>::const_iterator& parBgrOrderIt,
1688                      const int muonType) override {
1689     double thisStep[] = {0.01, 0.01, 0.01};
1690     TString thisParName[] = {"Bgr fraction", "Constant", "Linear"};
1691     if (muonType == 1) {
1692       double thisMini[] = {0.0, 0., -300.};
1693       double thisMaxi[] = {1.0, 300., 0.};
1694       this->setPar(
1695           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1696     } else {
1697       double thisMini[] = {0.0, 0., -300.};
1698       double thisMaxi[] = {1.0, 300., 0.};
1699       this->setPar(
1700           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1701     }
1702   }
1703 };
1704 
1705 /// Exponential binned in eta
1706 // --------------------------
1707 class backgroundFunctionType6 : public backgroundFunctionBase {
1708 public:
1709   backgroundFunctionType6(const double& lowerLimit, const double& upperLimit)
1710       : backgroundFunctionBase(lowerLimit, upperLimit) {
1711     this->parNum_ = 2;
1712   }
1713   double operator()(const double* parval, const double& mass, const double& eta) const override { return 0.; }
1714   double operator()(const double* parval, const double& mass, const double& eta1, const double& eta2) const override {
1715     double Bgrp2 = 0.;
1716     if (fabs(eta1) <= 1.3 && fabs(eta2) <= 1.3) {
1717       Bgrp2 = -1.20528;
1718     } else if ((fabs(eta1) <= 1.6 && fabs(eta1) > 1.3) && (fabs(eta2) <= 1.6 && fabs(eta2) > 1.3)) {
1719       Bgrp2 = 0.234713;
1720     } else if (fabs(eta1) > 1.6 && fabs(eta2) > 1.6) {
1721       Bgrp2 = -0.667103;
1722     } else if ((fabs(eta1) <= 1.3 && (fabs(eta2) > 1.3 && fabs(eta2) <= 1.6)) ||
1723                (fabs(eta2) <= 1.3 && (fabs(eta1) > 1.3 && fabs(eta1) <= 1.6))) {
1724       Bgrp2 = -0.656904;
1725     } else if ((fabs(eta1) <= 1.3 && fabs(eta2) > 1.6) || (fabs(eta2) <= 1.3 && fabs(eta1) > 1.6)) {
1726       Bgrp2 = 0.155328;
1727     } else if (((fabs(eta1) > 1.3 && fabs(eta1) <= 1.6) && fabs(eta2) > 1.6) ||
1728                ((fabs(eta2) > 1.3 && fabs(eta2) <= 1.6) && fabs(eta1) > 1.6)) {
1729       Bgrp2 = -0.177154;
1730     } else {
1731       std::cout << "WARNING: this should not happen for eta1 = " << eta1 << " and eta2 = " << eta2 << std::endl;
1732       Bgrp2 = -0.667103;
1733     }
1734     Bgrp2 *= -1.;
1735     double norm = -(exp(-Bgrp2 * upperLimit_) - exp(-Bgrp2 * lowerLimit_)) / Bgrp2;
1736     if (norm != 0)
1737       return exp(-Bgrp2 * mass) / norm;
1738     else
1739       return 0.;
1740   }
1741   void setParameters(double* Start,
1742                      double* Step,
1743                      double* Mini,
1744                      double* Maxi,
1745                      int* ind,
1746                      TString* parname,
1747                      const std::vector<double>::const_iterator& parBgrIt,
1748                      const std::vector<int>::const_iterator& parBgrOrderIt,
1749                      const int muonType) override {
1750     double thisStep[] = {0.01, 0.01};
1751     TString thisParName[] = {"Bgr fraction", "Bgr slope"};
1752     if (muonType == 1) {
1753       double thisMini[] = {0.0, 0.};
1754       double thisMaxi[] = {1.0, 10.};
1755       this->setPar(
1756           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1757     } else {
1758       double thisMini[] = {0.0, 0.};
1759       double thisMaxi[] = {1.0, 10.};
1760       this->setPar(
1761           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1762     }
1763   }
1764 
1765   double fracVsEta(const double* parval, const double& eta1, const double& eta2) const override {
1766     if (fabs(eta1) <= 1.3 && fabs(eta2) <= 1.3) {
1767       return (1. - 0.910903);
1768     } else if ((fabs(eta1) <= 1.6 && fabs(eta1) > 1.3) && (fabs(eta2) <= 1.6 && fabs(eta2) > 1.3)) {
1769       return (1. - 0.801469);
1770     } else if (fabs(eta1) > 1.6 && fabs(eta2) > 1.6) {
1771       return (1. - 0.658196);
1772     } else if ((fabs(eta1) <= 1.3 && (fabs(eta2) > 1.3 && fabs(eta2) <= 1.6)) ||
1773                (fabs(eta2) <= 1.3 && (fabs(eta1) > 1.3 && fabs(eta1) <= 1.6))) {
1774       return (1. - 0.873411);
1775     } else if ((fabs(eta1) <= 1.3 && fabs(eta2) > 1.6) || (fabs(eta2) <= 1.3 && fabs(eta1) > 1.6)) {
1776       return (1. - 0.784674);
1777     } else if (((fabs(eta1) > 1.3 && fabs(eta1) <= 1.6) && fabs(eta2) > 1.6) ||
1778                ((fabs(eta2) > 1.3 && fabs(eta2) <= 1.6) && fabs(eta1) > 1.6)) {
1779       return (1. - 0.714398);
1780     } else {
1781       std::cout << "WARNING: this should not happen for eta1 = " << eta1 << " and eta2 = " << eta2 << std::endl;
1782       return (1. - 0.658196);
1783     }
1784   }
1785 };
1786 
1787 /// Exponential binned in eta, much finer binning then type6
1788 // ---------------------------------------------------------
1789 class backgroundFunctionType7 : public backgroundFunctionBase {
1790 public:
1791   backgroundFunctionType7(const double& lowerLimit, const double& upperLimit)
1792       : backgroundFunctionBase(lowerLimit, upperLimit) {
1793     this->parNum_ = 2;
1794   }
1795   double operator()(const double* parval, const double& mass, const double& eta) const override { return 0.; }
1796   double operator()(const double* parval, const double& mass, const double& eta1, const double& eta2) const override {
1797     double Bgrp2 = 0.;
1798     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 0. && fabs(eta2) < 0.9)) {
1799       Bgrp2 = (-1.42465);
1800     } else if ((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 0.9 && fabs(eta2) < 1.3)) {
1801       Bgrp2 = (-1.38576);
1802     } else if ((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.3 && fabs(eta2) < 1.5)) {
1803       Bgrp2 = (-0.333728);
1804     } else if ((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) {
1805       Bgrp2 = (0.94066);
1806     } else if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) {
1807       Bgrp2 = (0.371026);
1808     } else if ((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) {
1809       Bgrp2 = (-0.959101);
1810     } else if ((fabs(eta1) >= 1.8 && fabs(eta1) < 1.9) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) {
1811       Bgrp2 = (-1.13829);
1812     } else if ((fabs(eta1) >= 1.9 && fabs(eta1) < 2.0) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) {
1813       Bgrp2 = (-0.921581);
1814     } else if ((fabs(eta1) >= 2.0 && fabs(eta1) < 1000.) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) {
1815       Bgrp2 = (-0.664338);
1816     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 0.9 && fabs(eta2) < 1.3)) ||
1817                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 0.9 && fabs(eta1) < 1.3))) {
1818       Bgrp2 = (-1.07581);
1819     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.3 && fabs(eta2) < 1.5)) ||
1820                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.3 && fabs(eta1) < 1.5))) {
1821       Bgrp2 = (-0.250272);
1822     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) ||
1823                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.5 && fabs(eta1) < 1.6))) {
1824       Bgrp2 = (0.101785);
1825     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
1826                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
1827       Bgrp2 = (0.360397);
1828     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
1829                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
1830       Bgrp2 = (0.689136);
1831     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
1832                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
1833       Bgrp2 = (0.860723);
1834     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
1835                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
1836       Bgrp2 = (1.21908);
1837     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1838                ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1839       Bgrp2 = (2.4453);
1840     } else if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.3 && fabs(eta2) < 1.5)) ||
1841                ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.3 && fabs(eta1) < 1.5))) {
1842       Bgrp2 = (-1.14152);
1843     } else if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) ||
1844                ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.5 && fabs(eta1) < 1.6))) {
1845       Bgrp2 = (-0.77241);
1846     } else if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
1847                ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
1848       Bgrp2 = (-0.516479);
1849     } else if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
1850                ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
1851       Bgrp2 = (-0.361401);
1852     } else if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
1853                ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
1854       Bgrp2 = (-0.33143);
1855     } else if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
1856                ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
1857       Bgrp2 = (-0.20813);
1858     } else if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1859                ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1860       Bgrp2 = (0.158002);
1861     } else if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) ||
1862                ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.5 && fabs(eta1) < 1.6))) {
1863       Bgrp2 = (0.273222);
1864     } else if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
1865                ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
1866       Bgrp2 = (0.247639);
1867     } else if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
1868                ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
1869       Bgrp2 = (-0.148616);
1870     } else if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
1871                ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
1872       Bgrp2 = (-0.413175);
1873     } else if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
1874                ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
1875       Bgrp2 = (-0.230031);
1876     } else if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1877                ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1878       Bgrp2 = (-0.122756);
1879     } else if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
1880                ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
1881       Bgrp2 = (0.650851);
1882     } else if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
1883                ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
1884       Bgrp2 = (-0.0985001);
1885     } else if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
1886                ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
1887       Bgrp2 = (-0.402548);
1888     } else if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
1889                ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
1890       Bgrp2 = (-0.27401);
1891     } else if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1892                ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1893       Bgrp2 = (-0.22863);
1894     } else if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
1895                ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
1896       Bgrp2 = (-0.436959);
1897     } else if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
1898                ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
1899       Bgrp2 = (-0.506041);
1900     } else if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
1901                ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
1902       Bgrp2 = (-0.31618);
1903     } else if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1904                ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1905       Bgrp2 = (-0.365653);
1906     } else if (((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
1907                ((fabs(eta2) >= 1.7 && fabs(eta2) < 1.8) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
1908       Bgrp2 = (-1.16783);
1909     } else if (((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
1910                ((fabs(eta2) >= 1.7 && fabs(eta2) < 1.8) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
1911       Bgrp2 = (-0.730701);
1912     } else if (((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1913                ((fabs(eta2) >= 1.7 && fabs(eta2) < 1.8) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1914       Bgrp2 = (-0.5271);
1915     } else if (((fabs(eta1) >= 1.8 && fabs(eta1) < 1.9) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
1916                ((fabs(eta2) >= 1.8 && fabs(eta2) < 1.9) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
1917       Bgrp2 = (-0.99893);
1918     } else if (((fabs(eta1) >= 1.8 && fabs(eta1) < 1.9) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1919                ((fabs(eta2) >= 1.8 && fabs(eta2) < 1.9) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1920       Bgrp2 = (-0.687263);
1921     } else if (((fabs(eta1) >= 1.9 && fabs(eta1) < 2.0) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
1922                ((fabs(eta2) >= 1.9 && fabs(eta2) < 2.0) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
1923       Bgrp2 = (-0.722394);
1924     } else {
1925       std::cout << "WARNING: this should not happen for eta1 = " << eta1 << " and eta2 = " << eta2 << std::endl;
1926       Bgrp2 = -0.664338;
1927     }
1928     Bgrp2 *= -1.;
1929     double norm = -(exp(-Bgrp2 * upperLimit_) - exp(-Bgrp2 * lowerLimit_)) / Bgrp2;
1930     if (norm != 0)
1931       return exp(-Bgrp2 * mass) / norm;
1932     else
1933       return 0.;
1934   }
1935   void setParameters(double* Start,
1936                      double* Step,
1937                      double* Mini,
1938                      double* Maxi,
1939                      int* ind,
1940                      TString* parname,
1941                      const std::vector<double>::const_iterator& parBgrIt,
1942                      const std::vector<int>::const_iterator& parBgrOrderIt,
1943                      const int muonType) override {
1944     double thisStep[] = {0.01, 0.01};
1945     TString thisParName[] = {"Bgr fraction", "Bgr slope"};
1946     if (muonType == 1) {
1947       double thisMini[] = {0.0, 0.};
1948       double thisMaxi[] = {1.0, 10.};
1949       this->setPar(
1950           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1951     } else {
1952       double thisMini[] = {0.0, 0.};
1953       double thisMaxi[] = {1.0, 10.};
1954       this->setPar(
1955           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
1956     }
1957   }
1958 
1959   double fracVsEta(const double* parval, const double& eta1, const double& eta2) const override {
1960     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 0. && fabs(eta2) < 0.9)) {
1961       return (1. - 0.915365);
1962     }
1963     if ((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 0.9 && fabs(eta2) < 1.3)) {
1964       return (1. - 0.914149);
1965     }
1966     if ((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.3 && fabs(eta2) < 1.5)) {
1967       return (1. - 0.855918);
1968     }
1969     if ((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) {
1970       return (1. - 0.70221);
1971     }
1972     if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) {
1973       return (1. - 0.701489);
1974     }
1975     if ((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) {
1976       return (1. - 0.651162);
1977     }
1978     if ((fabs(eta1) >= 1.8 && fabs(eta1) < 1.9) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) {
1979       return (1. - 0.639839);
1980     }
1981     if ((fabs(eta1) >= 1.9 && fabs(eta1) < 2.0) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) {
1982       return (1. - 0.64915);
1983     }
1984     if ((fabs(eta1) >= 2.0 && fabs(eta1) < 1000.) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) {
1985       return (1. - 0.687878);
1986     }
1987     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 0.9 && fabs(eta2) < 1.3)) ||
1988         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 0.9 && fabs(eta1) < 1.3))) {
1989       return (1. - 0.903486);
1990     }
1991     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.3 && fabs(eta2) < 1.5)) ||
1992         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.3 && fabs(eta1) < 1.5))) {
1993       return (1. - 0.882516);
1994     }
1995     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) ||
1996         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.5 && fabs(eta1) < 1.6))) {
1997       return (1. - 0.85477);
1998     }
1999     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
2000         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
2001       return (1. - 0.804919);
2002     }
2003     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
2004         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
2005       return (1. - 0.75411);
2006     }
2007     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
2008         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
2009       return (1. - 0.714128);
2010     }
2011     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
2012         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
2013       return (1. - 0.645403);
2014     }
2015     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.9) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2016         ((fabs(eta2) >= 0. && fabs(eta2) < 0.9) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2017       return (1. - 0.588049);
2018     }
2019     if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.3 && fabs(eta2) < 1.5)) ||
2020         ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.3 && fabs(eta1) < 1.5))) {
2021       return (1. - 0.901123);
2022     }
2023     if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) ||
2024         ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.5 && fabs(eta1) < 1.6))) {
2025       return (1. - 0.87852);
2026     }
2027     if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
2028         ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
2029       return (1. - 0.862266);
2030     }
2031     if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
2032         ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
2033       return (1. - 0.846385);
2034     }
2035     if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
2036         ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
2037       return (1. - 0.825401);
2038     }
2039     if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
2040         ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
2041       return (1. - 0.812449);
2042     }
2043     if (((fabs(eta1) >= 0.9 && fabs(eta1) < 1.3) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2044         ((fabs(eta2) >= 0.9 && fabs(eta2) < 1.3) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2045       return (1. - 0.753754);
2046     }
2047     if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.5 && fabs(eta2) < 1.6)) ||
2048         ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.5 && fabs(eta1) < 1.6))) {
2049       return (1. - 0.794143);
2050     }
2051     if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
2052         ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
2053       return (1. - 0.761375);
2054     }
2055     if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
2056         ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
2057       return (1. - 0.765572);
2058     }
2059     if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
2060         ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
2061       return (1. - 0.749438);
2062     }
2063     if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
2064         ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
2065       return (1. - 0.750941);
2066     }
2067     if (((fabs(eta1) >= 1.3 && fabs(eta1) < 1.5) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2068         ((fabs(eta2) >= 1.3 && fabs(eta2) < 1.5) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2069       return (1. - 0.722832);
2070     }
2071     if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.7)) ||
2072         ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.7))) {
2073       return (1. - 0.699723);
2074     }
2075     if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
2076         ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
2077       return (1. - 0.734044);
2078     }
2079     if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
2080         ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
2081       return (1. - 0.719434);
2082     }
2083     if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
2084         ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
2085       return (1. - 0.718889);
2086     }
2087     if (((fabs(eta1) >= 1.5 && fabs(eta1) < 1.6) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2088         ((fabs(eta2) >= 1.5 && fabs(eta2) < 1.6) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2089       return (1. - 0.689382);
2090     }
2091     if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.7 && fabs(eta2) < 1.8)) ||
2092         ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 1.7 && fabs(eta1) < 1.8))) {
2093       return (1. - 0.681106);
2094     }
2095     if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
2096         ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
2097       return (1. - 0.685783);
2098     }
2099     if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
2100         ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
2101       return (1. - 0.695924);
2102     }
2103     if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.7) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2104         ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.7) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2105       return (1. - 0.670977);
2106     }
2107     if (((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.8 && fabs(eta2) < 1.9)) ||
2108         ((fabs(eta2) >= 1.7 && fabs(eta2) < 1.8) && (fabs(eta1) >= 1.8 && fabs(eta1) < 1.9))) {
2109       return (1. - 0.654816);
2110     }
2111     if (((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
2112         ((fabs(eta2) >= 1.7 && fabs(eta2) < 1.8) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
2113       return (1. - 0.670969);
2114     }
2115     if (((fabs(eta1) >= 1.7 && fabs(eta1) < 1.8) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2116         ((fabs(eta2) >= 1.7 && fabs(eta2) < 1.8) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2117       return (1. - 0.659082);
2118     }
2119     if (((fabs(eta1) >= 1.8 && fabs(eta1) < 1.9) && (fabs(eta2) >= 1.9 && fabs(eta2) < 2.0)) ||
2120         ((fabs(eta2) >= 1.8 && fabs(eta2) < 1.9) && (fabs(eta1) >= 1.9 && fabs(eta1) < 2.0))) {
2121       return (1. - 0.648371);
2122     }
2123     if (((fabs(eta1) >= 1.8 && fabs(eta1) < 1.9) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2124         ((fabs(eta2) >= 1.8 && fabs(eta2) < 1.9) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2125       return (1. - 0.659114);
2126     }
2127     if (((fabs(eta1) >= 1.9 && fabs(eta1) < 2.0) && (fabs(eta2) >= 2.0 && fabs(eta2) < 1000.)) ||
2128         ((fabs(eta2) >= 1.9 && fabs(eta2) < 2.0) && (fabs(eta1) >= 2.0 && fabs(eta1) < 1000.))) {
2129       return (1. - 0.660482);
2130     } else {
2131       std::cout << "WARNING: this should not happen for eta1 = " << eta1 << " and eta2 = " << eta2 << std::endl;
2132       return (1. - 0.687878);
2133     }
2134   }
2135 };
2136 //Function 8
2137 class backgroundFunctionType8 : public backgroundFunctionBase {
2138 public:
2139   backgroundFunctionType8(const double& lowerLimit, const double& upperLimit)
2140       : backgroundFunctionBase(lowerLimit, upperLimit) {
2141     this->parNum_ = 2;
2142   }
2143   double operator()(const double* parval, const double& mass, const double& eta) const override { return 0.; }
2144   double operator()(const double* parval, const double& mass, const double& eta1, const double& eta2) const override {
2145     double Bgrp2 = 0.;
2146     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0. && fabs(eta2) < 0.85)) {
2147       Bgrp2 = (-2.17047);
2148     } else if ((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) {
2149       Bgrp2 = (-2.12913);
2150     } else if ((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) {
2151       Bgrp2 = (-2.19963);
2152     } else if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1000) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) {
2153       Bgrp2 = (-0.386394);
2154     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) ||
2155                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 0.85 && fabs(eta1) < 1.25))) {
2156       Bgrp2 = (-1.71339);
2157     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2158                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2159       Bgrp2 = (-0.206566);
2160     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2161                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2162       Bgrp2 = (4.4815);
2163     } else if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2164                ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2165       Bgrp2 = (-1.87985);
2166     } else if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2167                ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2168       Bgrp2 = (-0.163569);
2169     } else if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2170                ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2171       Bgrp2 = (-1.67935);
2172     }
2173     Bgrp2 *= -1.;
2174     double norm = -(exp(-Bgrp2 * upperLimit_) - exp(-Bgrp2 * lowerLimit_)) / Bgrp2;
2175     if (norm != 0)
2176       return exp(-Bgrp2 * mass) / norm;
2177     else
2178       return 0.;
2179   }
2180   void setParameters(double* Start,
2181                      double* Step,
2182                      double* Mini,
2183                      double* Maxi,
2184                      int* ind,
2185                      TString* parname,
2186                      const std::vector<double>::const_iterator& parBgrIt,
2187                      const std::vector<int>::const_iterator& parBgrOrderIt,
2188                      const int muonType) override {
2189     double thisStep[] = {0.01, 0.01};
2190     TString thisParName[] = {"Bgr fraction", "Bgr slope"};
2191     if (muonType == 1) {
2192       double thisMini[] = {0.0, 0.};
2193       double thisMaxi[] = {1.0, 10.};
2194       this->setPar(
2195           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2196     } else {
2197       double thisMini[] = {0.0, 0.};
2198       double thisMaxi[] = {1.0, 10.};
2199       this->setPar(
2200           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2201     }
2202   }
2203 
2204   double fracVsEta(const double* parval, const double& eta1, const double& eta2) const override {
2205     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0. && fabs(eta2) < 0.85)) {
2206       return (1. - 0.907727);
2207     }
2208     if ((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) {
2209       return (1. - 0.907715);
2210     }
2211     if ((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) {
2212       return (1. - 0.912233);
2213     }
2214     if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1000) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) {
2215       return (1. - 0.876776);
2216     }
2217     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) ||
2218         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 0.85 && fabs(eta1) < 1.25))) {
2219       return (1. - 0.913046);
2220     }
2221     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2222         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2223       return (1. - 0.916765);
2224     }
2225     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2226         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2227       return (1. - 0.6);
2228     }
2229     if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2230         ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2231       return (1. - 0.907471);
2232     }
2233     if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2234         ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2235       return (1. - 0.899253);
2236     }
2237     if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2238         ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2239       return (1. - 0.879108);
2240     }
2241     return 0;
2242   }
2243 };
2244 
2245 /// Service to build the background functor corresponding to the passed identifier
2246 backgroundFunctionBase* backgroundFunctionService(const int identifier,
2247                                                   const double& lowerLimit,
2248                                                   const double& upperLimit);
2249 
2250 //Function Type 9
2251 
2252 class backgroundFunctionType9 : public backgroundFunctionBase {
2253 public:
2254   backgroundFunctionType9(const double& lowerLimit, const double& upperLimit)
2255       : backgroundFunctionBase(lowerLimit, upperLimit) {
2256     this->parNum_ = 2;
2257   }
2258   double operator()(const double* parval, const double& mass, const double& eta) const override { return 0.; }
2259   double operator()(const double* parval, const double& mass, const double& eta1, const double& eta2) const override {
2260     double Bgrp2 = 0.;
2261 
2262     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0. && fabs(eta2) < 0.85)) {
2263       Bgrp2 = (-1.80833);
2264     } else if ((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) {
2265       Bgrp2 = (-1.98281);
2266     } else if ((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) {
2267       Bgrp2 = (-1.79632);
2268     } else if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1000) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) {
2269       Bgrp2 = (-1.14645);
2270     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) ||
2271                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 0.85 && fabs(eta1) < 1.25))) {
2272       Bgrp2 = (-1.55747);
2273     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2274                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2275       Bgrp2 = (-0.337598);
2276     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2277                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2278       Bgrp2 = (5.36513);
2279     } else if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2280                ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2281       Bgrp2 = (-1.44363);
2282     } else if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2283                ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2284       Bgrp2 = (-0.54614);
2285     } else if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2286                ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2287       Bgrp2 = (-1.41059);
2288     } else if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2289                ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2290       Bgrp2 = (-1.41059);
2291     }
2292     Bgrp2 *= -1.;
2293     double norm = -(exp(-Bgrp2 * upperLimit_) - exp(-Bgrp2 * lowerLimit_)) / Bgrp2;
2294     if (norm != 0)
2295       return exp(-Bgrp2 * mass) / norm;
2296     else
2297       return 0.;
2298   }
2299   void setParameters(double* Start,
2300                      double* Step,
2301                      double* Mini,
2302                      double* Maxi,
2303                      int* ind,
2304                      TString* parname,
2305                      const std::vector<double>::const_iterator& parBgrIt,
2306                      const std::vector<int>::const_iterator& parBgrOrderIt,
2307                      const int muonType) override {
2308     double thisStep[] = {0.01, 0.01};
2309     TString thisParName[] = {"Bgr fraction", "Bgr slope"};
2310     if (muonType == 1) {
2311       double thisMini[] = {0.0, 0.};
2312       double thisMaxi[] = {1.0, 10.};
2313       this->setPar(
2314           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2315     } else {
2316       double thisMini[] = {0.0, 0.};
2317       double thisMaxi[] = {1.0, 10.};
2318       this->setPar(
2319           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2320     }
2321   }
2322 
2323   double fracVsEta(const double* parval, const double& eta1, const double& eta2) const override
2324 
2325   {
2326     // std::cout << "Eta1 " << eta1 << " eta2 = " << eta2 << std::endl;
2327     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0. && fabs(eta2) < 0.85)) {
2328       return (1. - 0.893683);
2329     }
2330     if ((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) {
2331       return (1. - 0.888968);
2332     }
2333     if ((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) {
2334       return (1. - 0.885926);
2335     }
2336     if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1000) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) {
2337       return (1. - 0.866615);
2338     }
2339     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) ||
2340         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 0.85 && fabs(eta1) < 1.25))) {
2341       return (1. - 0.892856);
2342     }
2343     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2344         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2345       return (1. - 0.884864);
2346     }
2347     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2348         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2349       return (1. - 0.6);
2350     }
2351     if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2352         ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2353       return (1. - 0.894739);
2354     }
2355     if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2356         ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2357       return (1. - 0.880597);
2358     }
2359     if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2360         ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2361       return (1. - 0.869165);
2362     }
2363     if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2364         ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2365       return (1. - 0.869165);
2366     }
2367     return 0;
2368   }
2369 };
2370 
2371 class backgroundFunctionType10 : public backgroundFunctionBase {
2372 public:
2373   backgroundFunctionType10(const double& lowerLimit, const double& upperLimit)
2374       : backgroundFunctionBase(lowerLimit, upperLimit) {
2375     this->parNum_ = 2;
2376   }
2377   double operator()(const double* parval, const double& mass, const double& eta) const override { return 0.; }
2378   double operator()(const double* parval, const double& mass, const double& eta1, const double& eta2) const override {
2379     double Bgrp2 = 0.;
2380     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0. && fabs(eta2) < 0.85)) {
2381       Bgrp2 = (-1.80833);
2382     } else if ((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) {
2383       Bgrp2 = (-1.98281);
2384     } else if ((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) {
2385       Bgrp2 = (-1.79632);
2386     } else if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.8)) {
2387       Bgrp2 = (-1.50855);
2388     } else if ((fabs(eta1) >= 1.8 && fabs(eta1) < 2.0) && (fabs(eta2) >= 1.8 && fabs(eta2) < 2.0)) {
2389       Bgrp2 = (-0.498511);
2390     } else if ((fabs(eta1) >= 2.0 && fabs(eta1) < 2.2) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) {
2391       Bgrp2 = (-0.897031);
2392     } else if ((fabs(eta1) >= 2.2 && fabs(eta1) < 1000) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) {
2393       Bgrp2 = (-0.75954);
2394     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) ||
2395                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 0.85 && fabs(eta1) < 1.25))) {
2396       Bgrp2 = (-1.55747);
2397     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2398                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2399       Bgrp2 = (-0.337598);
2400     } else if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2401                ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2402       Bgrp2 = (3.5163);
2403     } else if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2404                ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2405       Bgrp2 = (-1.44363);
2406     } else if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2407                ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2408       Bgrp2 = (-0.54614);
2409     } else if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.8)) ||
2410                ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.8))) {
2411       Bgrp2 = (-1.36442);
2412     } else if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.8 && fabs(eta2) < 2.0)) ||
2413                ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.8 && fabs(eta1) < 2.0))) {
2414       Bgrp2 = (-1.66202);
2415     } else if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) ||
2416                ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 2.0 && fabs(eta1) < 2.2))) {
2417       Bgrp2 = (-0.62038);
2418     } else if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2419                ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2420       Bgrp2 = (0.662449);
2421     } else if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.8 && fabs(eta2) < 2.0)) ||
2422                ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.8) && (fabs(eta1) >= 1.8 && fabs(eta1) < 2.0))) {
2423       Bgrp2 = (-0.723325);
2424     } else if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) ||
2425                ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.8) && (fabs(eta1) >= 2.0 && fabs(eta1) < 2.2))) {
2426       Bgrp2 = (-1.54405);
2427     } else if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2428                ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.8) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2429       Bgrp2 = (-1.1104);
2430     } else if (((fabs(eta1) >= 1.8 && fabs(eta1) < 2.0) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) ||
2431                ((fabs(eta2) >= 1.8 && fabs(eta2) < 2.0) && (fabs(eta1) >= 2.0 && fabs(eta1) < 2.2))) {
2432       Bgrp2 = (-1.56277);
2433     } else if (((fabs(eta1) >= 1.8 && fabs(eta1) < 2.0) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2434                ((fabs(eta2) >= 1.8 && fabs(eta2) < 2.0) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2435       Bgrp2 = (-1.0827);
2436     } else if (((fabs(eta1) >= 2.0 && fabs(eta1) < 2.2) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2437                ((fabs(eta2) >= 2.0 && fabs(eta2) < 2.2) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2438       Bgrp2 = (-1.05662);
2439     }
2440     Bgrp2 *= -1.;
2441     double norm = -(exp(-Bgrp2 * upperLimit_) - exp(-Bgrp2 * lowerLimit_)) / Bgrp2;
2442     if (norm != 0)
2443       return exp(-Bgrp2 * mass) / norm;
2444     else
2445       return 0.;
2446   }
2447   void setParameters(double* Start,
2448                      double* Step,
2449                      double* Mini,
2450                      double* Maxi,
2451                      int* ind,
2452                      TString* parname,
2453                      const std::vector<double>::const_iterator& parBgrIt,
2454                      const std::vector<int>::const_iterator& parBgrOrderIt,
2455                      const int muonType) override {
2456     double thisStep[] = {0.01, 0.01};
2457     TString thisParName[] = {"Bgr fraction", "Bgr slope"};
2458     if (muonType == 1) {
2459       double thisMini[] = {0.0, 0.};
2460       double thisMaxi[] = {1.0, 10.};
2461       this->setPar(
2462           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2463     } else {
2464       double thisMini[] = {0.0, 0.};
2465       double thisMaxi[] = {1.0, 10.};
2466       this->setPar(
2467           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2468     }
2469   }
2470 
2471   double fracVsEta(const double* parval, const double& eta1, const double& eta2) const override {
2472     if ((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0. && fabs(eta2) < 0.85)) {
2473       return (1. - 0.893683);
2474     }
2475     if ((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) {
2476       return (1. - 0.888968);
2477     }
2478     if ((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) {
2479       return (1. - 0.885926);
2480     }
2481     if ((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.8)) {
2482       return (1. - 0.892823);
2483     }
2484     if ((fabs(eta1) >= 1.8 && fabs(eta1) < 2.0) && (fabs(eta2) >= 1.8 && fabs(eta2) < 2.0)) {
2485       return (1. - 0.888735);
2486     }
2487     if ((fabs(eta1) >= 2.0 && fabs(eta1) < 2.2) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) {
2488       return (1. - 0.87497);
2489     }
2490     if ((fabs(eta1) >= 2.2 && fabs(eta1) < 1000) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) {
2491       return (1. - 0.895275);
2492     }
2493     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 0.85 && fabs(eta2) < 1.25)) ||
2494         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 0.85 && fabs(eta1) < 1.25))) {
2495       return (1. - 0.892856);
2496     }
2497     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2498         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2499       return (1. - 0.884864);
2500     }
2501     if (((fabs(eta1) >= 0. && fabs(eta1) < 0.85) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2502         ((fabs(eta2) >= 0. && fabs(eta2) < 0.85) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2503       return (1. - 0.834572);
2504     }
2505     if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.25 && fabs(eta2) < 1.6)) ||
2506         ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.25 && fabs(eta1) < 1.6))) {
2507       return (1. - 0.894739);
2508     }
2509     if (((fabs(eta1) >= 0.85 && fabs(eta1) < 1.25) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1000)) ||
2510         ((fabs(eta2) >= 0.85 && fabs(eta2) < 1.25) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1000))) {
2511       return (1. - 0.880597);
2512     }
2513     if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.6 && fabs(eta2) < 1.8)) ||
2514         ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.6 && fabs(eta1) < 1.8))) {
2515       return (1. - 0.892911);
2516     }
2517     if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 1.8 && fabs(eta2) < 2.0)) ||
2518         ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 1.8 && fabs(eta1) < 2.0))) {
2519       return (1. - 0.880506);
2520     }
2521     if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) ||
2522         ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 2.0 && fabs(eta1) < 2.2))) {
2523       return (1. - 0.885718);
2524     }
2525     if (((fabs(eta1) >= 1.25 && fabs(eta1) < 1.6) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2526         ((fabs(eta2) >= 1.25 && fabs(eta2) < 1.6) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2527       return (1. - 0.853141);
2528     }
2529     if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 1.8 && fabs(eta2) < 2.0)) ||
2530         ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.8) && (fabs(eta1) >= 1.8 && fabs(eta1) < 2.0))) {
2531       return (1. - 0.88822);
2532     }
2533     if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) ||
2534         ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.8) && (fabs(eta1) >= 2.0 && fabs(eta1) < 2.2))) {
2535       return (1. - 0.87028);
2536     }
2537     if (((fabs(eta1) >= 1.6 && fabs(eta1) < 1.8) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2538         ((fabs(eta2) >= 1.6 && fabs(eta2) < 1.8) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2539       return (1. - 0.869603);
2540     }
2541     if (((fabs(eta1) >= 1.8 && fabs(eta1) < 2.0) && (fabs(eta2) >= 2.0 && fabs(eta2) < 2.2)) ||
2542         ((fabs(eta2) >= 1.8 && fabs(eta2) < 2.0) && (fabs(eta1) >= 2.0 && fabs(eta1) < 2.2))) {
2543       return (1. - 0.877922);
2544     }
2545     if (((fabs(eta1) >= 1.8 && fabs(eta1) < 2.0) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2546         ((fabs(eta2) >= 1.8 && fabs(eta2) < 2.0) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2547       return (1. - 0.865997);
2548     }
2549     if (((fabs(eta1) >= 2.0 && fabs(eta1) < 2.2) && (fabs(eta2) >= 2.2 && fabs(eta2) < 1000)) ||
2550         ((fabs(eta2) >= 2.0 && fabs(eta2) < 2.2) && (fabs(eta1) >= 2.2 && fabs(eta1) < 1000))) {
2551       return (1. - 0.886109);
2552     }
2553     return 0;
2554   }
2555 };
2556 
2557 /// Exponential binned in eta (Z, Run2012C PromptReco-v1 + PromptReco-v2)
2558 // --------------------------
2559 class backgroundFunctionType11 : public backgroundFunctionBase {
2560 public:
2561   backgroundFunctionType11(const double& lowerLimit, const double& upperLimit)
2562       : backgroundFunctionBase(lowerLimit, upperLimit) {
2563     this->parNum_ = 2;
2564   }
2565   double operator()(const double* parval, const double& mass, const double& eta) const override { return 0.; }
2566   double operator()(const double* parval, const double& mass, const double& eta1, const double& eta2) const override {
2567     double Bgrp2 = 0.;
2568     if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= -100. && eta2 < -0.8)) {
2569       Bgrp2 = (-0.0512353);
2570     } else if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= -0.8 && eta2 < 0.)) {
2571       Bgrp2 = (-0.0448482);
2572     } else if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= 0. && eta2 < 0.8)) {
2573       Bgrp2 = (-0.0193726);
2574     } else if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= 0.8 && eta2 < 100.)) {
2575       Bgrp2 = (0.0225765);
2576     } else if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= -100. && eta2 < -0.8)) {
2577       Bgrp2 = (-0.0822936);
2578     } else if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= -0.8 && eta2 < 0.)) {
2579       Bgrp2 = (-0.0676357);
2580     } else if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= 0. && eta2 < 0.8)) {
2581       Bgrp2 = (-0.0591544);
2582     } else if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= 0.8 && eta2 < 100.)) {
2583       Bgrp2 = (-0.0235858);
2584     } else if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= -100. && eta2 < -0.8)) {
2585       Bgrp2 = (-0.0317051);
2586     } else if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= -0.8 && eta2 < 0.)) {
2587       Bgrp2 = (-0.06139);
2588     } else if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= 0. && eta2 < 0.8)) {
2589       Bgrp2 = (-0.0747737);
2590     } else if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= 0.8 && eta2 < 100.)) {
2591       Bgrp2 = (-0.0810139);
2592     } else if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= -100. && eta2 < -0.8)) {
2593       Bgrp2 = (0.0229602);
2594     } else if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= -0.8 && eta2 < 0.)) {
2595       Bgrp2 = (-0.0224212);
2596     } else if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= 0. && eta2 < 0.8)) {
2597       Bgrp2 = (-0.0446273);
2598     } else if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= 0.8 && eta2 < 100.)) {
2599       Bgrp2 = (-0.0554561);
2600     } else {
2601       std::cout << "WARNING, backgroundFunctionType11: this should not happen for eta1 = " << eta1
2602                 << " and eta2 = " << eta2 << std::endl;
2603       return (-0.05);
2604     }
2605     double norm = (exp(Bgrp2 * upperLimit_) - exp(Bgrp2 * lowerLimit_)) / Bgrp2;
2606     if (norm != 0)
2607       return exp(Bgrp2 * mass) / norm;
2608     else
2609       return 0.;
2610   }
2611   void setParameters(double* Start,
2612                      double* Step,
2613                      double* Mini,
2614                      double* Maxi,
2615                      int* ind,
2616                      TString* parname,
2617                      const std::vector<double>::const_iterator& parBgrIt,
2618                      const std::vector<int>::const_iterator& parBgrOrderIt,
2619                      const int muonType) override {
2620     double thisStep[] = {0.01, 0.01};
2621     TString thisParName[] = {"Bgr fraction", "Bgr slope"};
2622     if (muonType == 1) {
2623       double thisMini[] = {-1.0, 10.};
2624       double thisMaxi[] = {1.0, 10.};
2625       this->setPar(
2626           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2627     } else {
2628       double thisMini[] = {-1.0, 10.};
2629       double thisMaxi[] = {1.0, 10.};
2630       this->setPar(
2631           Start, Step, Mini, Maxi, ind, parname, parBgrIt, parBgrOrderIt, thisStep, thisMini, thisMaxi, thisParName);
2632     }
2633   }
2634 
2635   double fracVsEta(const double* parval, const double& eta1, const double& eta2) const override {
2636     if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= -100. && eta2 < -0.8)) {
2637       return (1. - 0.966316);
2638     }
2639     if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= -0.8 && eta2 < 0.)) {
2640       return (1. - 0.966875);
2641     }
2642     if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= 0. && eta2 < 0.8)) {
2643       return (1. - 0.955311);
2644     }
2645     if ((eta1 >= -100. && eta1 < -0.8) && (eta2 >= 0.8 && eta2 < 100.)) {
2646       return (1. - 0.928771);
2647     }
2648     if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= -100. && eta2 < -0.8)) {
2649       return (1. - 0.983255);
2650     }
2651     if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= -0.8 && eta2 < 0.)) {
2652       return (1. - 0.982203);
2653     }
2654     if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= 0. && eta2 < 0.8)) {
2655       return (1. - 0.972127);
2656     }
2657     if ((eta1 >= -0.8 && eta1 < 0.) && (eta2 >= 0.8 && eta2 < 100.)) {
2658       return (1. - 0.962929);
2659     }
2660     if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= -100. && eta2 < -0.8)) {
2661       return (1. - 0.965597);
2662     }
2663     if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= -0.8 && eta2 < 0.)) {
2664       return (1. - 0.969461);
2665     }
2666     if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= 0. && eta2 < 0.8)) {
2667       return (1. - 0.979922);
2668     }
2669     if ((eta1 >= 0. && eta1 < 0.8) && (eta2 >= 0.8 && eta2 < 100.)) {
2670       return (1. - 0.984247);
2671     }
2672     if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= -100. && eta2 < -0.8)) {
2673       return (1. - 0.934252);
2674     }
2675     if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= -0.8 && eta2 < 0.)) {
2676       return (1. - 0.952914);
2677     }
2678     if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= 0. && eta2 < 0.8)) {
2679       return (1. - 0.960191);
2680     }
2681     if ((eta1 >= 0.8 && eta1 < 100.) && (eta2 >= 0.8 && eta2 < 100.)) {
2682       return (1. - 0.966175);
2683     } else {
2684       std::cout << "WARNING, backgroundFunctionType11: this should not happen for eta1 = " << eta1
2685                 << " and eta2 = " << eta2 << std::endl;
2686       return (1. - 0.97);
2687     }
2688   }
2689 };
2690 
2691 #endif  // FUNCTIONS_H