File indexing completed on 2024-04-06 12:26:08
0001
0002
0003
0004
0005
0006
0007 #include "RecoLocalMuon/DTSegment/src/DTLinearFit.h"
0008
0009
0010
0011
0012 #include <cmath>
0013 using namespace std;
0014
0015
0016
0017
0018 DTLinearFit::DTLinearFit() {}
0019
0020
0021 DTLinearFit::~DTLinearFit() {}
0022
0023
0024 void DTLinearFit::fit(const vector<float>& x,
0025 const vector<float>& y,
0026 int nptfit,
0027 const vector<float>& sigy,
0028 float& slope,
0029 float& intercept,
0030 double& chi2,
0031 float& covss,
0032 float& covii,
0033 float& covsi) const {
0034 chi2 = 0;
0035 float sy = 0, sx = 0;
0036 float s11 = 0, sxy = 0, sxx = 0;
0037 for (int i = 0; i != nptfit; i++) {
0038 float sy2 = sigy[i] * sigy[i];
0039 sy += y[i] / sy2;
0040 sxy += x[i] * y[i] / sy2;
0041 s11 += 1. / sy2;
0042 sx += x[i] / sy2;
0043 sxx += x[i] * x[i] / sy2;
0044 }
0045
0046 float delta = s11 * sxx - sx * sx;
0047
0048 if (delta == 0)
0049 return;
0050
0051 intercept = (sy * sxx - sxy * sx) / delta;
0052 slope = (sxy * s11 - sy * sx) / delta;
0053
0054 covii = sxx / delta;
0055 covss = s11 / delta;
0056 covsi = -sx / delta;
0057
0058 for (int j = 0; j < nptfit; j++) {
0059 const double ypred = intercept + slope * x[j];
0060 const double dy = (y[j] - ypred) / sigy[j];
0061 chi2 += dy * dy;
0062 }
0063 }
0064
0065 void DTLinearFit::fitNpar(const int npar,
0066 const vector<float>& xfit,
0067 const vector<float>& yfit,
0068 const vector<int>& lfit,
0069 const vector<double>& tfit,
0070 const vector<float>& sigy,
0071 float& aminf,
0072 float& bminf,
0073 float& cminf,
0074 float& vminf,
0075 double& chi2fit,
0076 const bool debug = false) const {
0077 int nptfit = xfit.size();
0078 if (nptfit < npar)
0079 return;
0080
0081 double sx = 0.;
0082 double sxx = 0.;
0083 double sy = 0.;
0084 double sxy = 0.;
0085 double sl = 0.;
0086 double sl2 = 0.;
0087 double sly = 0.;
0088 double slx = 0.;
0089 double st = 0.;
0090 double st2 = 0.;
0091 double slt = 0.;
0092 double sltx = 0.;
0093 double slty = 0.;
0094
0095 for (int j = 0; j < nptfit; j++) {
0096 sx += xfit[j];
0097 sy += yfit[j];
0098 sxx += xfit[j] * xfit[j];
0099 sxy += xfit[j] * yfit[j];
0100 sl += lfit[j];
0101 sl2 += lfit[j] * lfit[j];
0102 sly += lfit[j] * yfit[j];
0103 slx += lfit[j] * xfit[j];
0104 st += tfit[j];
0105 st2 += tfit[j] * tfit[j];
0106 slt += lfit[j] * tfit[j];
0107 sltx += lfit[j] * tfit[j] * xfit[j];
0108 slty += lfit[j] * tfit[j] * yfit[j];
0109 }
0110
0111 const double d1 = sy;
0112 const double d2 = sxy;
0113 const double d3 = sly;
0114 const double c1 = sl;
0115 const double c2 = slx;
0116 const double c3 = sl2;
0117 const double b1 = sx;
0118 const double b2 = sxx;
0119 const double b3 = slx;
0120 const double a1 = nptfit;
0121 const double a2 = sx;
0122 const double a3 = sl;
0123
0124 const double b4 = b2 * a1 - b1 * a2;
0125 const double c4 = c2 * a1 - c1 * a2;
0126 const double d4 = d2 * a1 - d1 * a2;
0127 const double b5 = a1 * b3 - a3 * b1;
0128 const double c5 = a1 * c3 - a3 * c1;
0129 const double d5 = a1 * d3 - d1 * a3;
0130 const double a6 = slt;
0131 const double b6 = sltx;
0132 const double c6 = st;
0133 const double v6 = st2;
0134 const double d6 = slty;
0135
0136 chi2fit = -1.;
0137
0138 if (npar == 3) {
0139 if (((c5 * b4 - c4 * b5) * b4 * a1) != 0) {
0140 cminf = (d5 * b4 - d4 * b5) / (c5 * b4 - c4 * b5);
0141 if (fabs(cminf) < 0.000001)
0142 cminf = 0;
0143 aminf = d4 / b4 - cminf * c4 / b4;
0144 bminf = (d1 / a1 - cminf * c1 / a1 - aminf * b1 / a1);
0145
0146 chi2fit = 0.;
0147 for (int j = 0; j < nptfit; j++) {
0148 const double ypred = aminf * xfit[j] + bminf;
0149 const double ycorr = yfit[j] - cminf * lfit[j];
0150 const double dy = (ycorr - ypred) / sigy[j];
0151 chi2fit += dy * dy;
0152
0153 if (tfit[j] == 0.)
0154 continue;
0155 double xwire = yfit[j] - tfit[j];
0156
0157
0158 double margin = 0.1;
0159
0160
0161
0162 if ((fabs(yfit[j] - xwire) > margin) && (fabs(ycorr - xwire) > margin) &&
0163 ((yfit[j] - xwire) * (ycorr - xwire) < 0)) {
0164
0165
0166 chi2fit = -1.;
0167 return;
0168 }
0169
0170 if (fabs(ypred - xwire) > 2.1 + margin) {
0171
0172
0173 chi2fit = -1.;
0174 return;
0175 }
0176
0177 if (fabs(ycorr - xwire) > 2.1 + margin) {
0178
0179
0180 chi2fit = -1.;
0181 return;
0182 }
0183 if (fabs(chi2fit < 0.0001))
0184 chi2fit = 0;
0185 }
0186 }
0187 } else if (npar == 4) {
0188 const double det =
0189 (a1 * a1 * (b2 * v6 - b6 * b6) -
0190 a1 * (a2 * a2 * v6 - 2 * a2 * a6 * b6 + a6 * a6 * b2 + b2 * c6 * c6 + b3 * (b3 * v6 - 2 * b6 * c6)) +
0191 a2 * a2 * c6 * c6 + 2 * a2 * (a3 * (b3 * v6 - b6 * c6) - a6 * b3 * c6) + a3 * a3 * (b6 * b6 - b2 * v6) +
0192 a6 * (2 * a3 * (b2 * c6 - b3 * b6) + a6 * b3 * b3));
0193
0194 if (det != 0) {
0195
0196 bminf = (a1 * (a2 * (b6 * d6 - v6 * d2) + a6 * (b6 * d2 - b2 * d6) + d1 * (b2 * v6 - b6 * b6)) -
0197 a2 * (b3 * (c6 * d6 - v6 * d3) + c6 * (b6 * d3 - c6 * d2)) +
0198 a3 * (b2 * (c6 * d6 - v6 * d3) + b3 * (v6 * d2 - b6 * d6) + b6 * (b6 * d3 - c6 * d2)) +
0199 a6 * (b2 * c6 * d3 + b3 * (b3 * d6 - b6 * d3 - c6 * d2)) -
0200 d1 * (b2 * c6 * c6 + b3 * (b3 * v6 - 2 * b6 * c6))) /
0201 det;
0202 aminf =
0203 -(a1 * a1 * (b6 * d6 - v6 * d2) -
0204 a1 * (a2 * (a6 * d6 - v6 * d1) - a6 * a6 * d2 + a6 * b6 * d1 + b3 * (c6 * d6 - v6 * d3) +
0205 c6 * (b6 * d3 - c6 * d2)) +
0206 a2 * (a3 * (c6 * d6 - v6 * d3) + c6 * (a6 * d3 - c6 * d1)) + a3 * a3 * (v6 * d2 - b6 * d6) +
0207 a3 * (a6 * (b3 * d6 + b6 * d3 - 2 * c6 * d2) - d1 * (b3 * v6 - b6 * c6)) - a6 * b3 * (a6 * d3 - c6 * d1)) /
0208 det;
0209 cminf = -(a1 * (b2 * (c6 * d6 - v6 * d3) + b3 * (v6 * d2 - b6 * d6) + b6 * (b6 * d3 - c6 * d2)) +
0210 a2 * a2 * (v6 * d3 - c6 * d6) +
0211 a2 * (a3 * (b6 * d6 - v6 * d2) + a6 * (b3 * d6 - 2 * b6 * d3 + c6 * d2) - d1 * (b3 * v6 - b6 * c6)) +
0212 a3 * (d1 * (b2 * v6 - b6 * b6) - a6 * (b2 * d6 - b6 * d2)) +
0213 a6 * (a6 * (b2 * d3 - b3 * d2) - d1 * (b2 * c6 - b3 * b6))) /
0214 det;
0215 vminf = -(a1 * a1 * (b2 * d6 - b6 * d2) -
0216 a1 * (a2 * a2 * d6 - a2 * (a6 * d2 + b6 * d1) + a6 * b2 * d1 + b2 * c6 * d3 +
0217 b3 * (b3 * d6 - b6 * d3 - c6 * d2)) +
0218 a2 * a2 * c6 * d3 + a2 * (a3 * (2 * b3 * d6 - b6 * d3 - c6 * d2) - b3 * (a6 * d3 + c6 * d1)) +
0219 a3 * a3 * (b6 * d2 - b2 * d6) + a3 * (a6 * (b2 * d3 - b3 * d2) + d1 * (b2 * c6 - b3 * b6)) +
0220 a6 * b3 * b3 * d1) /
0221 det;
0222
0223 chi2fit = 0.;
0224 for (int j = 0; j < nptfit; j++) {
0225 const double ypred = aminf * xfit[j] + bminf;
0226 const double dy = (yfit[j] + vminf * lfit[j] * tfit[j] - cminf * lfit[j] - ypred) / sigy[j];
0227 chi2fit += dy * dy;
0228 }
0229 }
0230 }
0231 }
0232
0233 void DTLinearFit::fit3par(const vector<float>& xfit,
0234 const vector<float>& yfit,
0235 const vector<int>& lfit,
0236 const int nptfit,
0237 const vector<float>& sigy,
0238 float& aminf,
0239 float& bminf,
0240 float& cminf,
0241 double& chi2fit,
0242 const bool debug = false) const {
0243 float vminf;
0244 vector<double> tfit(nptfit, 0.);
0245
0246 fitNpar(3, xfit, yfit, lfit, tfit, sigy, aminf, bminf, cminf, vminf, chi2fit, debug);
0247 }
0248
0249 void DTLinearFit::fit4Var(const vector<float>& xfit,
0250 const vector<float>& yfit,
0251 const vector<int>& lfit,
0252 const vector<double>& tfit,
0253 const int nptfit,
0254 float& aminf,
0255 float& bminf,
0256 float& cminf,
0257 float& vminf,
0258 double& chi2fit,
0259 const bool vdrift_4parfit = false,
0260 const bool debug = false) const {
0261 if (debug)
0262 cout << "Entering Fit4Var" << endl;
0263
0264 const double sigma =
0265 0.0295;
0266 vector<float> sigy;
0267
0268 aminf = 0.;
0269 bminf = 0.;
0270 cminf = -999.;
0271 vminf = 0.;
0272
0273 int nppar = 0;
0274 double chi2fitN2 = -1.;
0275 double chi2fit3 = -1.;
0276 double chi2fitN3 = -1.;
0277 double chi2fitN4 = -1.;
0278 int nppar2 = 0;
0279 int nppar3 = 0;
0280 int nppar4 = 0;
0281
0282 sigy.reserve(nptfit);
0283 for (int j = 0; j < nptfit; j++)
0284 sigy.push_back(sigma);
0285
0286 float a = 0.;
0287 float b = 0.;
0288 float covss, covii, covsi;
0289
0290 fit(xfit, yfit, nptfit, sigy, b, a, chi2fit, covss, covii, covsi);
0291
0292 bminf = b;
0293 aminf = a;
0294 nppar = 2;
0295 nppar2 = nppar;
0296 chi2fitN2 = chi2fit / (nptfit - nppar);
0297
0298
0299
0300 if (nptfit >= 3) {
0301 fitNpar(3, xfit, yfit, lfit, tfit, sigy, bminf, aminf, cminf, vminf, chi2fit, debug);
0302 chi2fit3 = chi2fit;
0303
0304 if (cminf != -999.) {
0305 nppar = 3;
0306 } else {
0307 bminf = b;
0308 aminf = a;
0309 }
0310 chi2fitN3 = (nptfit > nppar) ? chi2fit / (nptfit - nppar) : -1.;
0311
0312 float aminf3 = aminf;
0313 float bminf3 = bminf;
0314 float cminf3 = cminf;
0315 nppar3 = nppar;
0316
0317 if (debug) {
0318 cout << "dt0= 0 : slope 2 = " << b << " pos in = " << a << " chi2fitN2 = " << chi2fitN2 << " nppar = " << nppar2
0319 << " nptfit = " << nptfit << endl;
0320 cout << "dt0 = 0 : slope 3 = " << bminf << " pos out = " << aminf << " chi2fitN3 = " << chi2fitN3
0321 << " nppar = " << nppar3 << " T0_ev ns = " << cminf / 0.00543 << endl;
0322 cout << " vdrift_4parfit " << vdrift_4parfit << endl;
0323 }
0324
0325 if (nptfit >= 5) {
0326 fitNpar(4, xfit, yfit, lfit, tfit, sigy, bminf, aminf, cminf, vminf, chi2fit, debug);
0327
0328 if (vminf != 0)
0329 nppar = 4;
0330
0331 chi2fitN4 = (nptfit > nppar) ? chi2fit / (nptfit - nppar) : -1.;
0332
0333 if (fabs(vminf) >= 0.29) {
0334
0335 vminf = 0.;
0336 aminf = aminf3;
0337 bminf = bminf3;
0338 cminf = cminf3;
0339 nppar = 3;
0340 chi2fit = chi2fit3;
0341 }
0342 }
0343
0344 if (!vdrift_4parfit) {
0345
0346 aminf = aminf3;
0347 bminf = bminf3;
0348 cminf = cminf3;
0349 nppar = 3;
0350 chi2fit = chi2fit3;
0351 }
0352
0353 nppar4 = nppar;
0354
0355 }
0356
0357 if (debug) {
0358 cout << " dt0= 0 : slope 4 = " << bminf << " pos out = " << aminf << " chi2fitN4 = " << chi2fitN4
0359 << " nppar= " << nppar4 << " T0_ev ns= " << cminf / 0.00543 << " delta v = " << vminf << endl;
0360 cout << nptfit << " nptfit "
0361 << " end chi2fit = " << chi2fit / (nptfit - nppar) << " T0_ev ns= " << cminf / 0.00543
0362 << " delta v = " << vminf << endl;
0363 }
0364
0365 if (fabs(vminf) >= 0.09 && debug) {
0366
0367
0368
0369
0370
0371
0372
0373 cout << nptfit << " nptfit "
0374 << " end chi2fit = " << chi2fit << "T0_ev ns= " << cminf / 0.00543 << "delta v = " << vminf << endl;
0375 }
0376
0377 if (nptfit != nppar)
0378 chi2fit = chi2fit / (nptfit - nppar);
0379 }