Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:49

0001 /** \file
0002  *
0003  * \author Stefano Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>
0004  */
0005 
0006 /* This Class Header */
0007 #include "RecoLocalMuon/DTSegment/src/DTLinearFit.h"
0008 
0009 /* Collaborating Class Header */
0010 
0011 /* C++ Headers */
0012 #include <cmath>
0013 using namespace std;
0014 
0015 /* ====================================================================== */
0016 
0017 /// Constructor
0018 DTLinearFit::DTLinearFit() {}
0019 
0020 /// Destructor
0021 DTLinearFit::~DTLinearFit() {}
0022 
0023 /* Operations */
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   }  //end loop
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         // define a small safety margin
0158         double margin = 0.1;
0159         //          cout << " pred: " << ypred << "  hit: " << ycorr << "chi2: " << dy*dy << " t0: " << -cminf/0.00543 << endl;
0160 
0161         // Sanity checks - check that the hit after the fit is still withing the DT cell volume
0162         if ((fabs(yfit[j] - xwire) > margin) && (fabs(ycorr - xwire) > margin) &&
0163             ((yfit[j] - xwire) * (ycorr - xwire) < 0)) {
0164           //          cout << " segmentpred: " << ypred << "   hit: " << ycorr << "   xwire: " << xwire << "   yhit: " << yfit[j] << "   t0: " << -cminf/0.00543 << endl;
0165           //          cout << " XXX hit moved across wire!!!" << endl;
0166           chi2fit = -1.;
0167           return;
0168         }
0169 
0170         if (fabs(ypred - xwire) > 2.1 + margin) {
0171           //          cout << " segmentpred: " << ypred << "   hit: " << ycorr << "   xwire: " << xwire << "   yhit: " << yfit[j] << "   t0: " << -cminf/0.00543 << endl;
0172           //          cout << " XXX segment outside chamber!!!" << endl;
0173           chi2fit = -1.;
0174           return;
0175         }
0176 
0177         if (fabs(ycorr - xwire) > 2.1 + margin) {
0178           //          cout << " segmentpred: " << ypred << "   hit: " << ycorr << "   xwire: " << xwire << "   yhit: " << yfit[j] << "   t0: " << -cminf/0.00543 << endl;
0179           //          cout << " XXX hit outside chamber!!!" << endl;
0180           chi2fit = -1.;
0181           return;
0182         }
0183         if (fabs(chi2fit < 0.0001))
0184           chi2fit = 0;
0185       }  //end loop chi2
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       // computation of a, b, c and v
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       }  //end loop chi2
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;  // FIXME errors can be inserted .just load them/that is the usual TB resolution value for DT chambers
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   // cout << "dt0 = 0  chi2fit = " << chi2fit << "  slope = "<<b<<endl;
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         // for safety and for code construction..dont accept correction on dv/vdrift greater then 0.09
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) {  //if not required explicitly leave the t0 and track step as at step 3
0345                             // just update vdrift value vmin for storing in the segments for monitoring
0346       aminf = aminf3;
0347       bminf = bminf3;
0348       cminf = cminf3;
0349       nppar = 3;
0350       chi2fit = chi2fit3;
0351     }
0352 
0353     nppar4 = nppar;
0354 
0355   }  //end nptfit >=3
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) {  //checks only vdrift less then 10 % accepted
0366                                        //    cout << "vminf gt 0.09 det=  " << endl;
0367     //    cout << "dt0= 0 : slope 4 = "<< bminf << " pos out = " << aminf << " chi2fitN4 = " << chi2fitN4
0368     //   << " T0_ev ns = " << cminf/0.00543 << " delta v = "<< vminf << endl;
0369     //    cout << "dt0 = 0 : slope 2 = "<< b << " pos in = " << a <<" chi2fitN2 = " << chi2fitN2
0370     //   << " nppar = " << nppar-1 << " nptfit = " << nptfit <<endl;
0371     //    cout << "dt0 = 0 : slope 3 = " << bminf << " pos out = " << aminf << " chi2fitN3 = "
0372     //   << chi2fitN3 << " T0_ev ns = " << cminf/0.00543 << endl;
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 }