CSCSegment/src/CSCSegFit.cc

0001 // CSCSegFit.cc
0002 // Last mod: 03.02.2015
0003
0004 #include "RecoLocalMuon/CSCSegment/src/CSCSegFit.h"
0005
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007
0008 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 void CSCSegFit::fit(void) {
0014   if (fitdone())
0015     return;  // don't redo fit unnecessarily
0016   short n = nhits();
0017   switch (n) {
0018     case 1:
0019       edm::LogVerbatim("CSCSegFit") << "[CSCSegFit::fit] - cannot fit just 1 hit!!";
0020       break;
0021     case 2:
0022       fit2();
0023       break;
0024     case 3:
0025     case 4:
0026     case 5:
0027     case 6:
0028       fitlsq();
0029       break;
0030     default:
0031       edm::LogVerbatim("CSCSegFit") << "[CSCSegFit::fit] - cannot fit more than 6 hits!!";
0032   }
0033 }
0034
0035 void CSCSegFit::fit2(void) {
0036   // Just join the two points
0037   // Equation of straight line between (x1, y1) and (x2, y2) in xy-plane is
0038   //       y = mx + c
0039   // with m = (y2-y1)/(x2-x1)
0040   // and  c = (y1*x2-x2*y1)/(x2-x1)
0041
0042   CSCSetOfHits::const_iterator ih = hits_.begin();
0043   int il1 = (*ih)->cscDetId().layer();
0044   const CSCRecHit2D& h1 = (**ih);
0045   ++ih;
0046   int il2 = (*ih)->cscDetId().layer();
0047   const CSCRecHit2D& h2 = (**ih);
0048
0049   // Skip if on same layer, but should be impossible :)
0050   if (il1 == il2) {
0051     edm::LogVerbatim("CSCSegFit") << "[CSCSegFit:fit]2 - 2 hits on same layer!!";
0052     return;
0053   }
0054
0055   const CSCLayer* layer1 = chamber()->layer(il1);
0056   const CSCLayer* layer2 = chamber()->layer(il2);
0057
0058   GlobalPoint h1glopos = layer1->toGlobal(h1.localPosition());
0059   GlobalPoint h2glopos = layer2->toGlobal(h2.localPosition());
0060
0061   // We want hit wrt chamber (and local z will be != 0)
0062   LocalPoint h1pos = chamber()->toLocal(h1glopos);
0063   LocalPoint h2pos = chamber()->toLocal(h2glopos);
0064
0065   float dz = h2pos.z() - h1pos.z();
0066
0067   uslope_ = (h2pos.x() - h1pos.x()) / dz;
0068   vslope_ = (h2pos.y() - h1pos.y()) / dz;
0069
0070   float uintercept = (h1pos.x() * h2pos.z() - h2pos.x() * h1pos.z()) / dz;
0071   float vintercept = (h1pos.y() * h2pos.z() - h2pos.y() * h1pos.z()) / dz;
0072   intercept_ = LocalPoint(uintercept, vintercept, 0.);
0073
0074   setOutFromIP();
0075
0076   //@@ NOT SURE WHAT IS SENSIBLE FOR THESE...
0077   chi2_ = 0.;
0078   ndof_ = 0;
0079
0080   fitdone_ = true;
0081 }
0082
0083 void CSCSegFit::fitlsq(void) {
0084   // Linear least-squares fit to up to 6 CSC rechits, one per layer in a CSC.
0085   // Comments adapted from  mine in original  CSCSegAlgoSK algorithm.
0086
0087   // Fit to the local x, y rechit coordinates in z projection
0088   // The strip measurement controls the precision of x
0089   // The wire measurement controls the precision of y.
0090   // Typical precision: u (strip, sigma~200um), v (wire, sigma~1cm)
0091
0092   // Set up the normal equations for the least-squares fit as a matrix equation
0093
0094   // We have a vector of measurements m, which is a 2n x 1 dim matrix
0095   // The transpose mT is (u1, v1, u2, v2, ..., un, vn) where
0096   // ui is the strip-associated measurement and
0097   // vi is the wire-associated measurement
0098   // for a given rechit i.
0099
0100   // The fit is to
0101   // u = u0 + uz * z
0102   // v = v0 + vz * z
0103   // where u0, uz, v0, vz are the parameters to be obtained from the fit.
0104
0105   // These are contained in a vector p which is a 4x1 dim matrix, and
0106   // its transpose pT is (u0, v0, uz, vz). Note the ordering!
0107
0108   // The covariance matrix for each pair of measurements is 2 x 2 and
0109   // the inverse of this is the error matrix E.
0110   // The error matrix for the whole set of n measurements is a diagonal
0111   // matrix with diagonal elements the individual 2 x 2 error matrices
0112   // (because the inverse of a diagonal matrix is a diagonal matrix
0113   // with each element the inverse of the original.)
0114
0115   // In function 'weightMatrix()', the variable 'matrix' is filled with this
0116   // block-diagonal overall covariance matrix. Then 'matrix' is inverted to the
0117   // block-diagonal error matrix, and returned.
0118
0119   // Define the matrix A as
0120   //    1   0   z1  0
0121   //    0   1   0   z1
0122   //    1   0   z2  0
0123   //    0   1   0   z2
0124   //    ..  ..  ..  ..
0125   //    1   0   zn  0
0126   //    0   1   0   zn
0127
0128   // This matrix A is set up and returned by function 'derivativeMatrix()'.
0129
0130   // Then the normal equations are described by the matrix equation
0131   //
0132   //    (AT E A)p = (AT E)m
0133   //
0134   // where AT is the transpose of A.
0135
0136   // Call the combined matrix on the LHS, M, and that on the RHS, B:
0137   //     M p = B
0138
0139   // We solve this for the parameter vector, p.
0140   // The elements of M and B then involve sums over the hits
0141
0142   // The covariance matrix of the parameters is obtained by
0143   // (AT E A)^-1 calculated in 'covarianceMatrix()'.
0144
0145   // NOTE
0146   // We need local position of a RecHit w.r.t. the CHAMBER
0147   // and the RecHit itself only knows its local position w.r.t.
0148   // the LAYER, so we must explicitly transform global position.
0149
0150   SMatrix4 M;  // 4x4, init to 0
0151   SVector4 B;  // 4x1, init to 0;
0152
0153   CSCSetOfHits::const_iterator ih = hits_.begin();
0154
0155   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0156     const CSCRecHit2D& hit = (**ih);
0157     const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0158     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0159     LocalPoint lp = chamber()->toLocal(gp);
0160
0161     // Local position of hit w.r.t. chamber
0162     double u = lp.x();
0163     double v = lp.y();
0164     double z = lp.z();
0165
0166     // Covariance matrix of local errors
0167     SMatrixSym2 IC;  // 2x2, init to 0
0168
0169     IC(0, 0) = hit.localPositionError().xx();
0170     IC(1, 1) = hit.localPositionError().yy();
0171     //@@ NOT SURE WHICH OFF-DIAGONAL ELEMENT MUST BE DEFINED BUT (1,0) WORKS
0172     //@@ (and SMatrix enforces symmetry)
0173     IC(1, 0) = hit.localPositionError().xy();
0174     // IC(0,1) = IC(1,0);
0175
0176     // Invert covariance matrix (and trap if it fails!)
0177     bool ok = IC.Invert();
0178     if (!ok) {
0179       edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::fit] Failed to invert covariance matrix: \n" << IC;
0180       //      return ok;  //@@ SHOULD PASS THIS BACK TO CALLER?
0181     }
0182
0183     M(0, 0) += IC(0, 0);
0184     M(0, 1) += IC(0, 1);
0185     M(0, 2) += IC(0, 0) * z;
0186     M(0, 3) += IC(0, 1) * z;
0187     B(0) += u * IC(0, 0) + v * IC(0, 1);
0188
0189     M(1, 0) += IC(1, 0);
0190     M(1, 1) += IC(1, 1);
0191     M(1, 2) += IC(1, 0) * z;
0192     M(1, 3) += IC(1, 1) * z;
0193     B(1) += u * IC(1, 0) + v * IC(1, 1);
0194
0195     M(2, 0) += IC(0, 0) * z;
0196     M(2, 1) += IC(0, 1) * z;
0197     M(2, 2) += IC(0, 0) * z * z;
0198     M(2, 3) += IC(0, 1) * z * z;
0199     B(2) += (u * IC(0, 0) + v * IC(0, 1)) * z;
0200
0201     M(3, 0) += IC(1, 0) * z;
0202     M(3, 1) += IC(1, 1) * z;
0203     M(3, 2) += IC(1, 0) * z * z;
0204     M(3, 3) += IC(1, 1) * z * z;
0205     B(3) += (u * IC(1, 0) + v * IC(1, 1)) * z;
0206   }
0207
0208   SVector4 p;
0209   bool ok = M.Invert();
0210   if (!ok) {
0211     edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::fit] Failed to invert matrix: \n" << M;
0212     //    return ok; //@@ SHOULD PASS THIS BACK TO CALLER?
0213   } else {
0214     p = M * B;
0215   }
0216
0217   //  LogTrace("CSCSegFit") << "[CSCSegFit::fit] p = "
0218   //        << p(0) << ", " << p(1) << ", " << p(2) << ", " << p(3);
0219
0220   // fill member variables  (note origin has local z = 0)
0221   //  intercept_
0222   intercept_ = LocalPoint(p(0), p(1), 0.);
0223
0224   // localdir_ - set so segment points outwards from IP
0225   uslope_ = p(2);
0226   vslope_ = p(3);
0227   setOutFromIP();
0228
0229   // calculate chi2 of fit
0230   setChi2();
0231
0232   // flag fit has been done
0233   fitdone_ = true;
0234 }
0235
0236 void CSCSegFit::setChi2(void) {
0237   double chsq = 0.;
0238
0239   CSCSetOfHits::const_iterator ih;
0240   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0241     const CSCRecHit2D& hit = (**ih);
0242     const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0243     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0244     LocalPoint lp = chamber()->toLocal(gp);
0245
0246     double u = lp.x();
0247     double v = lp.y();
0248     double z = lp.z();
0249
0250     double du = intercept_.x() + uslope_ * z - u;
0251     double dv = intercept_.y() + vslope_ * z - v;
0252
0253     //    LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] u, v, z = " << u << ", " << v << ", " << z;
0254
0255     SMatrixSym2 IC;  // 2x2, init to 0
0256
0257     IC(0, 0) = hit.localPositionError().xx();
0258     //    IC(0,1) = hit.localPositionError().xy();
0259     IC(1, 0) = hit.localPositionError().xy();
0260     IC(1, 1) = hit.localPositionError().yy();
0261     //    IC(1,0) = IC(0,1);
0262
0263     //    LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC before = \n" << IC;
0264
0265     // Invert covariance matrix
0266     bool ok = IC.Invert();
0267     if (!ok) {
0268       edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::setChi2] Failed to invert covariance matrix: \n" << IC;
0269       //      return ok;
0270     }
0271     //    LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] IC after = \n" << IC;
0272     chsq += du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0273   }
0274
0275   // fill member variables
0276   chi2_ = chsq;
0277   ndof_ = 2. * hits_.size() - 4;
0278
0279   //  LogTrace("CSCSegFit") << "[CSCSegFit::setChi2] chi2 = " << chi2_ << "/" << ndof_ << " dof";
0280 }
0281
0282 CSCSegFit::SMatrixSym12 CSCSegFit::weightMatrix() {
0283   bool ok = true;
0284
0285   SMatrixSym12 matrix = ROOT::Math::SMatrixIdentity();  // 12x12, init to 1's on diag
0286
0287   int row = 0;
0288
0289   for (CSCSetOfHits::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
0290     const CSCRecHit2D& hit = (**it);
0291
0292     // Note scaleXError allows rescaling the x error if necessary
0293
0294     matrix(row, row) = scaleXError() * hit.localPositionError().xx();
0295     matrix(row, row + 1) = hit.localPositionError().xy();
0296     ++row;
0297     matrix(row, row - 1) = hit.localPositionError().xy();
0298     matrix(row, row) = hit.localPositionError().yy();
0299     ++row;
0300   }
0301
0302   ok = matrix.Invert();  // invert in place
0303   if (!ok) {
0304     edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::weightMatrix] Failed to invert matrix: \n" << matrix;
0305     //    return ok; //@@ SHOULD PASS THIS BACK TO CALLER?
0306   }
0307   return matrix;
0308 }
0309
0310 CSCSegFit::SMatrix12by4 CSCSegFit::derivativeMatrix() {
0311   SMatrix12by4 matrix;  // 12x4, init to 0
0312   int row = 0;
0313
0314   for (CSCSetOfHits::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
0315     const CSCRecHit2D& hit = (**it);
0316     const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0317     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0318     LocalPoint lp = chamber()->toLocal(gp);
0319     float z = lp.z();
0320
0321     matrix(row, 0) = 1.;
0322     matrix(row, 2) = z;
0323     ++row;
0324     matrix(row, 1) = 1.;
0325     matrix(row, 3) = z;
0326     ++row;
0327   }
0328   return matrix;
0329 }
0330
0331 void CSCSegFit::setOutFromIP() {
0332   // Set direction of segment to point from IP outwards
0333   // (Incorrect for particles not coming from IP, of course.)
0334
0335   double dxdz = uslope_;
0336   double dydz = vslope_;
0337   double dz = 1. / sqrt(1. + dxdz * dxdz + dydz * dydz);
0338   double dx = dz * dxdz;
0339   double dy = dz * dydz;
0340   LocalVector localDir(dx, dy, dz);
0341
0342   // localDir sometimes needs a sign flip
0343   // Examine its direction and origin in global z: to point outward
0344   // the localDir should always have same sign as global z...
0345
0346   double globalZpos = (chamber()->toGlobal(intercept_)).z();
0347   double globalZdir = (chamber()->toGlobal(localDir)).z();
0348   double directionSign = globalZpos * globalZdir;
0349   localdir_ = (directionSign * localDir).unit();
0350 }
0351
0352 AlgebraicSymMatrix CSCSegFit::covarianceMatrix() {
0353   SMatrixSym12 weights = weightMatrix();
0354   SMatrix12by4 A = derivativeMatrix();
0355   //  LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] weights matrix W: \n" << weights;
0356   //  LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] derivatives matrix A: \n" << A;
0357
0358   // (AT W A)^-1
0359   // e.g. See http://www.phys.ufl.edu/~avery/fitting.html, part I
0360
0361   bool ok;
0362   SMatrixSym4 result = ROOT::Math::SimilarityT(A, weights);
0363   //  LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] (AT W A): \n" << result;
0364   ok = result.Invert();  // inverts in place
0365   if (!ok) {
0366     edm::LogVerbatim("CSCSegment|CSCSegFit") << "[CSCSegFit::calculateError] Failed to invert matrix: \n" << result;
0367     //    return ok;  //@@ SHOULD PASS THIS BACK TO CALLER?
0368   }
0369   //  LogTrace("CSCSegFit") << "[CSCSegFit::covarianceMatrix] (AT W A)^-1: \n" << result;
0370
0371   // reorder components to match TrackingRecHit interface (CSCSegment isa TrackingRecHit)
0372   // i.e. slopes first, then positions
0373   AlgebraicSymMatrix flipped = flipErrors(result);
0374
0375   return flipped;
0376 }
0377
0378 AlgebraicSymMatrix CSCSegFit::flipErrors(const SMatrixSym4& a) {
0379   // The CSCSegment needs the error matrix re-arranged to match
0380   // parameters in order (uz, vz, u0, v0)
0381   // where uz, vz = slopes, u0, v0 = intercepts
0382
0383   //  LogTrace("CSCSegFit") << "[CSCSegFit::flipErrors] input: \n" << a;
0384
0385   AlgebraicSymMatrix hold(4, 0.);
0386
0387   for (short j = 0; j != 4; ++j) {
0388     for (short i = 0; i != 4; ++i) {
0389       hold(i + 1, j + 1) = a(i, j);  // SMatrix counts from 0, AlgebraicMatrix from 1
0390     }
0391   }
0392
0393   //  LogTrace("CSCSegFit") << "[CSCSegFit::flipErrors] after copy:";
0394   //  LogTrace("CSCSegFit") << "(" << hold(1,1) << "  " << hold(1,2) << "  " << hold(1,3) << "  " << hold(1,4);
0395   //  LogTrace("CSCSegFit") << " " << hold(2,1) << "  " << hold(2,2) << "  " << hold(2,3) << "  " << hold(2,4);
0396   //  LogTrace("CSCSegFit") << " " << hold(3,1) << "  " << hold(3,2) << "  " << hold(3,3) << "  " << hold(3,4);
0397   //  LogTrace("CSCSegFit") << " " << hold(4,1) << "  " << hold(4,2) << "  " << hold(4,3) << "  " << hold(4,4) << ")";
0398
0399   // errors on slopes into upper left
0400   hold(1, 1) = a(2, 2);
0401   hold(1, 2) = a(2, 3);
0402   hold(2, 1) = a(3, 2);
0403   hold(2, 2) = a(3, 3);
0404
0405   // errors on positions into lower right
0406   hold(3, 3) = a(0, 0);
0407   hold(3, 4) = a(0, 1);
0408   hold(4, 3) = a(1, 0);
0409   hold(4, 4) = a(1, 1);
0410
0411   // must also interchange off-diagonal elements of off-diagonal 2x2 submatrices
0412   hold(4, 1) = a(1, 2);
0413   hold(3, 2) = a(0, 3);
0414   hold(2, 3) = a(3, 0);  // = a(0,3)
0415   hold(1, 4) = a(2, 1);  // = a(1,2)
0416
0417   //  LogTrace("CSCSegFit") << "[CSCSegFit::flipErrors] after flip:";
0418   //  LogTrace("CSCSegFit") << "(" << hold(1,1) << "  " << hold(1,2) << "  " << hold(1,3) << "  " << hold(1,4);
0419   //  LogTrace("CSCSegFit") << " " << hold(2,1) << "  " << hold(2,2) << "  " << hold(2,3) << "  " << hold(2,4);
0420   //  LogTrace("CSCSegFit") << " " << hold(3,1) << "  " << hold(3,2) << "  " << hold(3,3) << "  " << hold(3,4);
0421   //  LogTrace("CSCSegFit") << " " << hold(4,1) << "  " << hold(4,2) << "  " << hold(4,3) << "  " << hold(4,4) << ")";
0422
0423   return hold;
0424 }
0425
0426 float CSCSegFit::xfit(float z) const {
0427   //@@ ADD THIS TO EACH ACCESSOR OF FIT RESULTS?
0428   //  if ( !fitdone() ) fit();
0429   return intercept_.x() + uslope_ * z;
0430 }
0431
0432 float CSCSegFit::yfit(float z) const { return intercept_.y() + vslope_ * z; }
0433
0434 float CSCSegFit::xdev(float x, float z) const { return intercept_.x() + uslope_ * z - x; }
0435
0436 float CSCSegFit::ydev(float y, float z) const { return intercept_.y() + vslope_ * z - y; }
0437
0438 float CSCSegFit::Rdev(float x, float y, float z) const {
0439   return sqrt(xdev(x, z) * xdev(x, z) + ydev(y, z) * ydev(y, z));
0440 }