Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:19:27 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!