MomentumScaleCalibration/interface/SigmaPtDiff.h

0001 #ifndef SigmaPtDiff_h
0002 #define SigmaPtDiff_h
0003
0004 #include <vector>
0005 #include <cmath>
0006
0007 class SigmaPt {
0008 public:
0009   SigmaPt(const std::vector<double>& parameters_, const std::vector<double>& errors_) {
0010     setParErr(parameters_, errors_);
0011   }
0012
0013   SigmaPt(){};
0014
0015   void setParErr(const std::vector<double>& parameters, const std::vector<double>& errors) {
0016     b_0 = parameters[0];
0017     b_1 = parameters[5];
0018     b_2 = parameters[1];
0019     b_3 = parameters[7];
0020     b_4 = parameters[8];
0021     sb_0 = errors[0];
0022     sb_1 = errors[5];
0023     sb_2 = errors[1];
0024     sb_3 = errors[7];
0025     sb_4 = errors[8];
0026     c = b_2 + b_3 * (b_0 - b_4) * (b_0 - b_4) - b_1 * b_0 * b_0;
0027   }
0028
0029   double operator()(const double& eta) {
0030     if (fabs(eta) <= b_0) {
0031       return (c + b_1 * eta * eta);
0032     }
0033     return (b_2 + b_3 * (fabs(eta) - b_4) * (fabs(eta) - b_4));
0034   }
0035   double sigma(const double& eta) {
0036     if (fabs(eta) <= b_0) {
0037       return sqrt((eta * eta - b_0 * b_0) * (eta * eta - b_0 * b_0) * sb_1 * sb_1 + sb_2 * sb_2 +
0038                   pow(b_0 - b_4, 4) * sb_3 * sb_3 + pow(-2 * b_3 * pow(b_0 - b_4, 2), 2) * sb_4 * sb_4 +
0039                   pow(2 * b_3 * (b_0 - b_4) - 2 * b_1 * b_0, 2) * sb_0 * sb_0);
0040     }
0041     return sqrt(sb_2 * sb_2 + pow(fabs(eta) - b_4, 4) * sb_3 * sb_3 +
0042                 pow(-2 * b_3 * pow(fabs(eta) - b_4, 2), 2) * sb_4 * sb_4);
0043   }
0044
0045 protected:
0046   double b_0;
0047   double b_1;
0048   double b_2;
0049   double b_3;
0050   double b_4;
0051   double c;
0052
0053   double sb_0;
0054   double sb_1;
0055   double sb_2;
0056   double sb_3;
0057   double sb_4;
0058 };
0059
0060 /// Returns ( sigmaPt/Pt(data)^2 - sigmaPt/Pt(MC)^2 )
0061 class SigmaPtDiff {
0062 public:
0063   SigmaPtDiff() {
0064     std::vector<double> parameters;
0065     std::vector<double> errors;
0066     parameters.push_back(1.66);
0067     parameters.push_back(0.021);
0068     parameters.push_back(0.);
0069     parameters.push_back(0.);
0070     parameters.push_back(0.);
0071     parameters.push_back(0.0058);
0072     parameters.push_back(0.);
0073     parameters.push_back(0.03);
0074     parameters.push_back(1.8);
0075     parameters.push_back(0.);
0076     parameters.push_back(0.);
0077     parameters.push_back(0.);
0078     parameters.push_back(0.);
0079     parameters.push_back(0.);
0080     parameters.push_back(0.);
0081     errors.push_back(0.09);
0082     errors.push_back(0.002);
0083     errors.push_back(0.);
0084     errors.push_back(0.);
0085     errors.push_back(0.);
0086     errors.push_back(0.0009);
0087     errors.push_back(0.);
0088     errors.push_back(0.03);
0089     errors.push_back(0.3);
0090     errors.push_back(0.);
0091     errors.push_back(0.);
0092     errors.push_back(0.);
0093     errors.push_back(0.);
0094     errors.push_back(0.);
0095     errors.push_back(0.);
0096
0097     sigmaPt.setParErr(parameters, errors);
0098   }
0099   double etaByPoints(const double& eta) {
0100     if (eta <= -2.2)
0101       return 0.0233989;
0102     else if (eta <= -2.0)
0103       return 0.0197057;
0104     else if (eta <= -1.8)
0105       return 0.014693;
0106     else if (eta <= -1.6)
0107       return 0.0146727;
0108     else if (eta <= -1.4)
0109       return 0.0141323;
0110     else if (eta <= -1.2)
0111       return 0.0159712;
0112     else if (eta <= -1.0)
0113       return 0.0117224;
0114     else if (eta <= -0.8)
0115       return 0.010726;
0116     else if (eta <= -0.6)
0117       return 0.0104777;
0118     else if (eta <= -0.4)
0119       return 0.00814458;
0120     else if (eta <= -0.2)
0121       return 0.00632501;
0122     else if (eta <= 0.0)
0123       return 0.00644172;
0124     else if (eta <= 0.2)
0125       return 0.00772645;
0126     else if (eta <= 0.4)
0127       return 0.010103;
0128     else if (eta <= 0.6)
0129       return 0.0099275;
0130     else if (eta <= 0.8)
0131       return 0.0100309;
0132     else if (eta <= 1.0)
0133       return 0.0125116;
0134     else if (eta <= 1.2)
0135       return 0.0147211;
0136     else if (eta <= 1.4)
0137       return 0.0151623;
0138     else if (eta <= 1.6)
0139       return 0.015259;
0140     else if (eta <= 1.8)
0141       return 0.014499;
0142     else if (eta <= 2.0)
0143       return 0.0165215;
0144     else if (eta <= 2.2)
0145       return 0.0212348;
0146     return 0.0227285;
0147   }
0148   // double squaredDiff(const double & eta, SigmaPt & sigmaPt)
0149   double squaredDiff(const double& eta) {
0150     double sigmaPtPlus = sigmaPt(eta) + sigmaPt.sigma(eta);
0151     double sigmaPtMinus = sigmaPt(eta) - sigmaPt.sigma(eta);
0152     if (fabs(sigmaPtPlus * sigmaPtPlus - etaByPoints(eta) * etaByPoints(eta)) >
0153         fabs(sigmaPtMinus * sigmaPtMinus - etaByPoints(eta) * etaByPoints(eta))) {
0154       return (fabs(sigmaPtPlus * sigmaPtPlus - etaByPoints(eta) * etaByPoints(eta)));
0155     }
0156     return (fabs(sigmaPtMinus * sigmaPtMinus - etaByPoints(eta) * etaByPoints(eta)));
0157   }
0158   SigmaPt sigmaPt;
0159 };
0160
0161 #endif