CSCSegment/src/CSCCondSegFit.cc

0001 // CSCCondSegFit.cc
0002 // Last mod: 23.01.2015
0003
0004 #include "RecoLocalMuon/CSCSegment/src/CSCCondSegFit.h"
0005
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include <algorithm>
0011
0012 void CSCCondSegFit::fit(bool condpass1, bool condpass2) {
0013   // See base class CSCSegFit::fit for detailed explanation of algorithm.
0014   // This is exactly the same, except it adds two conditioning passes
0015   // which can improve the robustness of the calculations.
0016
0017   lex_.clear();  // Vector of the error matrix (only xx)
0018
0019   if (condpass1 && !condpass2) {
0020     correctTheCovX();
0021     if (lex_.size() != hits_.size()) {
0022       LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] lex_.size()!=hits_.size() ALARM! Please inform CSC DPG "
0023                                  << std::endl;
0024     }
0025   }
0026
0027   SMatrix4 M;  // 4x4, init to 0
0028   SVector4 B;  // 4x1, init to 0;
0029
0030   CSCSetOfHits::const_iterator ih = hits_.begin();
0031
0032   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0033     const CSCRecHit2D& hit = (**ih);
0034     const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0035     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0036     LocalPoint lp = chamber()->toLocal(gp);
0037
0038     // Local position of hit w.r.t. chamber
0039     double u = lp.x();
0040     double v = lp.y();
0041     double z = lp.z();
0042
0043     // Covariance matrix of local errors
0044     SMatrixSym2 IC;  // 2x2, init to 0
0045
0046     // Adjust covariance matrix elements to improve numerical stability?
0047     if (condpass1 && !condpass2) {
0048       IC(0, 0) = lex_.at(ih - hits_.begin());
0049     } else {
0050       IC(0, 0) = hit.localPositionError().xx();
0051     }
0052     IC(1, 1) = hit.localPositionError().yy();
0053     //@@ NOT SURE WHICH OFF-DIAGONAL ELEMENT MUST BE DEFINED BUT (1,0) WORKS
0054     //@@ (and SMatrix enforces symmetry)
0055     IC(1, 0) = hit.localPositionError().xy();
0056     // IC(0,1) = IC(1,0);
0057
0058     // Correct the covariance matrix to improve numerical stability?
0059     if (condpass2) {
0060       correctTheCovMatrix(IC);
0061     }
0062
0063     // Invert covariance matrix (and trap if it fails!)
0064     bool ok = IC.Invert();
0065     if (!ok) {
0066       LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] Failed to invert covariance matrix: \n" << IC;
0067       //      return ok;  //@@ SHOULD PASS THIS BACK TO CALLER?
0068     }
0069
0070     M(0, 0) += IC(0, 0);
0071     M(0, 1) += IC(0, 1);
0072     M(0, 2) += IC(0, 0) * z;
0073     M(0, 3) += IC(0, 1) * z;
0074     B(0) += u * IC(0, 0) + v * IC(0, 1);
0075
0076     M(1, 0) += IC(1, 0);
0077     M(1, 1) += IC(1, 1);
0078     M(1, 2) += IC(1, 0) * z;
0079     M(1, 3) += IC(1, 1) * z;
0080     B(1) += u * IC(1, 0) + v * IC(1, 1);
0081
0082     M(2, 0) += IC(0, 0) * z;
0083     M(2, 1) += IC(0, 1) * z;
0084     M(2, 2) += IC(0, 0) * z * z;
0085     M(2, 3) += IC(0, 1) * z * z;
0086     B(2) += (u * IC(0, 0) + v * IC(0, 1)) * z;
0087
0088     M(3, 0) += IC(1, 0) * z;
0089     M(3, 1) += IC(1, 1) * z;
0090     M(3, 2) += IC(1, 0) * z * z;
0091     M(3, 3) += IC(1, 1) * z * z;
0092     B(3) += (u * IC(1, 0) + v * IC(1, 1)) * z;
0093   }
0094
0095   SVector4 p;
0096   bool ok = M.Invert();
0097   if (!ok) {
0098     LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] Failed to invert matrix: \n" << M;
0099     //    return ok; //@@ SHOULD PASS THIS BACK TO CALLER?
0100   } else {
0101     p = M * B;
0102   }
0103
0104   //  LogTrace("CSCSegFit") << "[CSCSegFit::fit] p = "
0105   //        << p(0) << ", " << p(1) << ", " << p(2) << ", " << p(3);
0106
0107   // fill member variables  (note origin has local z = 0)
0108   intercept_ = LocalPoint(p(0), p(1), 0.);
0109
0110   // localdir_
0111   uslope_ = p(2);
0112   vslope_ = p(3);
0113   setOutFromIP();
0114
0115   // calculate chi2 of fit
0116   setChi2(condpass1, condpass2);
0117 }
0118
0119 void CSCCondSegFit::setChi2(bool condpass1, bool condpass2) {
0120   double chsq = 0.;
0121
0122   CSCSetOfHits::const_iterator ih;
0123   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0124     const CSCRecHit2D& hit = (**ih);
0125     const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0126     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0127     LocalPoint lp = chamber()->toLocal(gp);
0128
0129     double u = lp.x();
0130     double v = lp.y();
0131     double z = lp.z();
0132
0133     double du = intercept_.x() + uslope_ * z - u;
0134     double dv = intercept_.y() + vslope_ * z - v;
0135
0136     //    LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] u, v, z = " << u << ", " << v << ", " << z;
0137
0138     SMatrixSym2 IC;  // 2x2, init to 0
0139
0140     if (condpass1 && !condpass2) {
0141       IC(0, 0) = lex_.at(ih - hits_.begin());
0142     } else {
0143       IC(0, 0) = hit.localPositionError().xx();
0144     }
0145     //    IC(0,1) = hit.localPositionError().xy();
0146     IC(1, 0) = hit.localPositionError().xy();
0147     IC(1, 1) = hit.localPositionError().yy();
0148     //    IC(1,0) = IC(0,1);
0149
0150     //    LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC before = \n" << IC;
0151
0152     /// Correct the cov matrix
0153     if (condpass2) {
0154       correctTheCovMatrix(IC);
0155     }
0156
0157     // Invert covariance matrix
0158     bool ok = IC.Invert();
0159     if (!ok) {
0160       LogTrace("CSCSegment|CSC") << "[CSCSegFit::setChi2] Failed to invert covariance matrix: \n" << IC;
0161       //      return ok;
0162     }
0163     //    LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC after = \n" << IC;
0164     chsq += du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0165   }
0166
0167   // fill member variables
0168   chi2_ = chsq;
0169   ndof_ = 2. * hits_.size() - 4;
0170
0171   //  LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] chi2 = " << chi2_ << "/" << ndof_ << " dof";
0172 }
0173
0174 void CSCCondSegFit::correctTheCovX(void) {
0175   std::vector<double> uu, vv, zz;  // Vectors of coordinates
0176
0177   lex_.clear();  // x component of local error for each hit
0178
0179   double sum_U_err = 0.0;
0180   double sum_Z_U_err = 0.0;
0181   double sum_Z2_U_err = 0.0;
0182   double sum_U_U_err = 0.0;
0183   double sum_UZ_U_err = 0.0;
0184   std::vector<double> chiUZind;
0185   std::vector<double>::iterator chiContribution;
0186   double chiUZ = 0.0;
0187   CSCSetOfHits::const_iterator ih = hits_.begin();
0188   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0189     const CSCRecHit2D& hit = (**ih);
0190     lex_.push_back(hit.localPositionError().xx());
0191
0192     const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0193     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0194     LocalPoint lp = chamber()->toLocal(gp);
0195     // Local position of hit w.r.t. chamber
0196     double u = lp.x();
0197     double v = lp.y();
0198     double z = lp.z();
0199     uu.push_back(u);
0200     vv.push_back(v);
0201     zz.push_back(z);
0202     // Prepare the sums for the standard linear fit
0203     sum_U_err += 1. / lex_.back();
0204     sum_Z_U_err += z / lex_.back();
0205     sum_Z2_U_err += (z * z) / lex_.back();
0206     sum_U_U_err += u / lex_.back();
0207     sum_UZ_U_err += (u * z) / lex_.back();
0208   }
0209
0210   /// Make a one dimensional fit in U-Z plane (i.e. chamber local x-z)
0211   /// U=U0+UZ*Z fit parameters
0212
0213   double denom = sum_U_err * sum_Z2_U_err - pow(sum_Z_U_err, 2);
0214   double U0 = (sum_Z2_U_err * sum_U_U_err - sum_Z_U_err * sum_UZ_U_err) / denom;
0215   double UZ = (sum_U_err * sum_UZ_U_err - sum_Z_U_err * sum_U_U_err) / denom;
0216
0217   /// Calculate the fit line trace
0218   /// Calculate one dimensional chi^2 and normalize the errors if needed
0219
0220   for (unsigned i = 0; i < uu.size(); ++i) {
0221     double uMean = U0 + UZ * zz[i];
0222     chiUZind.push_back((pow((uMean - uu[i]), 2)) / lex_[i]);
0223     chiUZ += (pow((uMean - uu[i]), 2)) / lex_[i];
0224   }
0225   chiUZ = chiUZ / (uu.size() - 2);
0226
0227   if (chiUZ >= chi2Norm_) {
0228     double chi2uCorrection = chiUZ / chi2Norm_;
0229     for (unsigned i = 0; i < uu.size(); ++i)
0230       lex_[i] = lex_[i] * chi2uCorrection;
0231     setScaleXError(chi2uCorrection);
0232   }
0233
0234   // Find most deviant hit
0235   chiContribution = max_element(chiUZind.begin(), chiUZind.end());
0236   worstHit_ = chiContribution - chiUZind.begin();
0237 }
0238
0239 void CSCCondSegFit::correctTheCovMatrix(SMatrixSym2& IC) {
0240   //double condNumberCorr1=0.0;
0241   double condNumberCorr2 = 0.0;
0242   double detCov = 0.0;
0243   double diag1 = 0.0;
0244   double diag2 = 0.0;
0245   double IC_01_corr = 0.0;
0246   double IC_00_corr = 0.0;
0247   if (!covToAnyNumberAll_) {
0248     //condNumberCorr1=condSeed1_*IC(1,1);
0249     condNumberCorr2 = condSeed2_ * IC(1, 1);
0250     diag1 = IC(0, 0) * IC(1, 1);
0251     diag2 = IC(0, 1) * IC(0, 1);
0252     detCov = fabs(diag1 - diag2);
0253     if ((diag1 < condNumberCorr2) && (diag2 < condNumberCorr2)) {
0254       if (covToAnyNumber_)
0255         IC(0, 1) = covAnyNumber_;
0256       else {
0257         IC_00_corr = condSeed1_ + fabs(IC(0, 1)) / IC(1, 1);
0258         IC(0, 0) = IC_00_corr;
0259       }
0260     }
0261
0262     if (((detCov < condNumberCorr2) && (diag1 > condNumberCorr2)) ||
0263         ((diag2 > condNumberCorr2) && (detCov < condNumberCorr2))) {
0264       if (covToAnyNumber_)
0265         IC(0, 1) = covAnyNumber_;
0266       else {
0267         IC_01_corr = sqrt(fabs(diag1 - condNumberCorr2));
0268         if (IC(0, 1) < 0)
0269           IC(0, 1) = (-1) * IC_01_corr;
0270         else
0271           IC(0, 1) = IC_01_corr;
0272       }
0273     }
0274   } else {
0275     IC(0, 1) = covAnyNumber_;
0276   }
0277   // With SMatrix formulation. the following might be unnecessary since it might
0278   // enforce the symmetry. But can't hurt to leave it (it WAS a real bug fix for
0279   // the original CLHEP formulation.)
0280   IC(1, 0) = IC(0, 1);  //@@ ADDED TO MAINTAIN SYMMETRY OF IC
0281 }