Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:11 

0001 // ------------------------- //
0002 // --> GEMCSCSegFit.cc
0003 // Created:  21.04.2015
0004 // --> Based on CSCSegFit.cc
0005 // with last mod: 03.02.2015
0006 // as found in 750pre2 rel
0007 // ------------------------- //
0008 
0009 #include "RecoLocalMuon/GEMCSCSegment/plugins/GEMCSCSegFit.h"
0010 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0011 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0012 
0013 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0014 #include "Geometry/GEMGeometry/interface/GEMEtaPartition.h"
0015 
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 
0019 void GEMCSCSegFit::fit(void) {
0020   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] - start the fitting fun :: nhits = " << nhits();
0021   if (fitdone())
0022     return;  // don't redo fit unnecessarily
0023   short n = nhits();
0024   switch (n) {
0025     case 1:
0026       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] - cannot fit just 1 hit!!";
0027       break;
0028     case 2:
0029       fit2();
0030       break;
0031     case 3:
0032     case 4:
0033     case 5:
0034     case 6:
0035     case 7:
0036     case 8:
0037       fitlsq();
0038       break;
0039     default:
0040       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] - cannot fit more than 8 hits!!";
0041   }
0042 }
0043 
0044 void GEMCSCSegFit::fit2(void) {
0045   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit2] - start fit2()";
0046   // Just join the two points
0047   // Equation of straight line between (x1, y1) and (x2, y2) in xy-plane is
0048   //       y = mx + c
0049   // with m = (y2-y1)/(x2-x1)
0050   // and  c = (y1*x2-x2*y1)/(x2-x1)
0051 
0052   // 1) Check whether hits are on the same layer
0053   // -------------------------------------------
0054   std::vector<const TrackingRecHit*>::const_iterator ih = hits_.begin();
0055   // layer numbering: GEM: (1,2) CSC (3,4,5,6,7,8)
0056   int il1 = 0, il2 = 0;
0057   DetId d1 = DetId((*ih)->rawId());
0058   // check whether first hit is GEM or CSC
0059   if (d1.subdetId() == MuonSubdetId::GEM) {
0060     il1 = GEMDetId(d1).layer();
0061   } else if (d1.subdetId() == MuonSubdetId::CSC) {
0062     il1 = CSCDetId(d1).layer() + 2;
0063   } else {
0064     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0065   }
0066   const TrackingRecHit& h1 = (**ih);
0067   ++ih;
0068   DetId d2 = DetId((*ih)->rawId());
0069   // check whether second hit is GEM or CSC
0070   if (d2.subdetId() == MuonSubdetId::GEM) {
0071     il2 = GEMDetId(d2).layer();
0072   } else if (d2.subdetId() == MuonSubdetId::CSC) {
0073     il2 = GEMDetId(d2).layer() + 2;
0074   } else {
0075     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0076   }
0077   const TrackingRecHit& h2 = (**ih);
0078   // Skip if on same layer, but should be impossible :)
0079   if (il1 == il2) {
0080     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - 2 hits on same layer!!";
0081     return;
0082   }
0083 
0084   // 2) Global Positions of hit 1 and 2 and
0085   //    Local  Positions of hit 1 and 2 w.r.t. reference CSC Chamber Frame
0086   // ---------------------------------------------------------------------
0087   GlobalPoint h1glopos, h2glopos;
0088   // global position hit 1
0089   if (d1.subdetId() == MuonSubdetId::GEM) {
0090     const GEMEtaPartition* roll1 = gemetapartition(GEMDetId(d1));
0091     h1glopos = roll1->toGlobal(h1.localPosition());
0092   } else if (d1.subdetId() == MuonSubdetId::CSC) {
0093     const CSCLayer* layer1 = cscchamber(CSCDetId(d1))->layer(CSCDetId(d1).layer());
0094     h1glopos = layer1->toGlobal(h1.localPosition());
0095   } else {
0096     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0097   }
0098   // global position hit 2
0099   if (d2.subdetId() == MuonSubdetId::GEM) {
0100     const GEMEtaPartition* roll2 = gemetapartition(GEMDetId(d2));
0101     h2glopos = roll2->toGlobal(h2.localPosition());
0102   } else if (d2.subdetId() == MuonSubdetId::CSC) {
0103     const CSCLayer* layer2 = cscchamber(CSCDetId(d2))->layer(CSCDetId(d2).layer());
0104     h2glopos = layer2->toGlobal(h2.localPosition());
0105   } else {
0106     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0107   }
0108   // local positions hit 1 and 2 w.r.t. ref CSC Chamber Frame
0109   // We want hit wrt chamber (and local z will be != 0)
0110   LocalPoint h1pos = refcscchamber()->toLocal(h1glopos);
0111   LocalPoint h2pos = refcscchamber()->toLocal(h2glopos);
0112 
0113   // 3) Now make straight line between the two points in local coords
0114   // ----------------------------------------------------------------
0115   float dz = h2pos.z() - h1pos.z();
0116   if (dz != 0) {
0117     uslope_ = (h2pos.x() - h1pos.x()) / dz;
0118     vslope_ = (h2pos.y() - h1pos.y()) / dz;
0119   }
0120   float uintercept = (h1pos.x() * h2pos.z() - h2pos.x() * h1pos.z()) / dz;
0121   float vintercept = (h1pos.y() * h2pos.z() - h2pos.y() * h1pos.z()) / dz;
0122   intercept_ = LocalPoint(uintercept, vintercept, 0.);
0123 
0124   setOutFromIP();
0125 
0126   //@@ NOT SURE WHAT IS SENSIBLE FOR THESE...
0127   chi2_ = 0.;
0128   ndof_ = 0;
0129 
0130   fitdone_ = true;
0131 }
0132 
0133 void GEMCSCSegFit::fitlsq(void) {
0134   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - start fitlsq";
0135   // Linear least-squares fit to up to 6 CSC rechits, one per layer in a CSC
0136   // and up to 2 GEM rechits in a GEM superchamber
0137   // (we can later later on go up to 4 hits in case of an overlapping chamber)
0138   // (... and if there is an overlap in the GEM chambers, than maybe we also
0139   //  have an overlap in the CSC chambers... maybe we can also benefit there)
0140 
0141   // Comments below taken from CSCSegFit algorithm.
0142   // Comments adapted from original  CSCSegAlgoSK algorithm.
0143 
0144   // Fit to the local x, y rechit coordinates in z projection
0145   // The CSC & GEM strip measurement controls the precision of x
0146   // The CSC wire measurement & GEM eta-partition controls the precision of y.
0147   // Typical precision CSC: u (strip, sigma~200um), v (wire, sigma~1cm)
0148   // Typical precision GEM: u (strip, sigma~250um), v (eta-part, sigma~10-20cm)
0149 
0150   // Set up the normal equations for the least-squares fit as a matrix equation
0151 
0152   // We have a vector of measurements m, which is a 2n x 1 dim matrix
0153   // The transpose mT is (u1, v1, u2, v2, ..., un, vn) where
0154   // ui is the strip-associated measurement and
0155   // vi is the wire-associated measurement
0156   // for a given rechit i.
0157 
0158   // The fit is to
0159   // u = u0 + uz * z
0160   // v = v0 + vz * z
0161   // where u0, uz, v0, vz are the parameters to be obtained from the fit.
0162 
0163   // These are contained in a vector p which is a 4x1 dim matrix, and
0164   // its transpose pT is (u0, v0, uz, vz). Note the ordering!
0165 
0166   // The covariance matrix for each pair of measurements is 2 x 2 and
0167   // the inverse of this is the error matrix E.
0168   // The error matrix for the whole set of n measurements is a diagonal
0169   // matrix with diagonal elements the individual 2 x 2 error matrices
0170   // (because the inverse of a diagonal matrix is a diagonal matrix
0171   // with each element the inverse of the original.)
0172 
0173   // In function 'weightMatrix()', the variable 'matrix' is filled with this
0174   // block-diagonal overall covariance matrix. Then 'matrix' is inverted to the
0175   // block-diagonal error matrix, and returned.
0176 
0177   // Define the matrix A as
0178   //    1   0   z1  0
0179   //    0   1   0   z1
0180   //    1   0   z2  0
0181   //    0   1   0   z2
0182   //    ..  ..  ..  ..
0183   //    1   0   zn  0
0184   //    0   1   0   zn
0185 
0186   // This matrix A is set up and returned by function 'derivativeMatrix()'.
0187 
0188   // Then the normal equations are described by the matrix equation
0189   //
0190   //    (AT E A)p = (AT E)m
0191   //
0192   // where AT is the transpose of A.
0193 
0194   // Call the combined matrix on the LHS, M, and that on the RHS, B:
0195   //     M p = B
0196 
0197   // We solve this for the parameter vector, p.
0198   // The elements of M and B then involve sums over the hits
0199 
0200   // The covariance matrix of the parameters is obtained by
0201   // (AT E A)^-1 calculated in 'covarianceMatrix()'.
0202 
0203   // NOTE
0204   // We need local position of a RecHit w.r.t. the CHAMBER
0205   // and the RecHit itself only knows its local position w.r.t.
0206   // the LAYER, so we must explicitly transform global position.
0207 
0208   SMatrix4 M;  // 4x4, init to 0
0209   SVector4 B;  // 4x1, init to 0;
0210 
0211   std::vector<const TrackingRecHit*>::const_iterator ih = hits_.begin();
0212 
0213   // LogDebug :: Loop over the TrackingRecHits and print the GEM and CSC Hits
0214   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0215     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - looping over TrackingRecHits";
0216     const TrackingRecHit& hit = (**ih);
0217     DetId d = DetId(hit.rawId());
0218     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit in detid (" << d.rawId() << ")";
0219     if (d.subdetId() == MuonSubdetId::GEM) {
0220       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit is a GEM Hit in detid ("
0221                                        << d.rawId() << ")";
0222       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - GEM DetId (" << GEMDetId(d.rawId()) << ")";
0223     } else if (d.subdetId() == MuonSubdetId::CSC) {
0224       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit is a CSC Hit in detid ("
0225                                        << d.rawId() << ")";
0226       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - CSC DetId (" << CSCDetId(d.rawId()) << ")";
0227     } else {
0228       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0229     }
0230   }
0231 
0232   // Loop over the TrackingRecHits and make small (2x2) matrices used to fill the blockdiagonal covariance matrix E^-1
0233   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0234     const TrackingRecHit& hit = (**ih);
0235     GlobalPoint gp;
0236     DetId d = DetId(hit.rawId());
0237     if (d.subdetId() == MuonSubdetId::GEM) {
0238       const GEMEtaPartition* roll = gemetapartition(GEMDetId(d));
0239       gp = roll->toGlobal(hit.localPosition());
0240     } else if (d.subdetId() == MuonSubdetId::CSC) {
0241       const CSCLayer* layer = cscchamber(CSCDetId(d))->layer(CSCDetId(d).layer());
0242       gp = layer->toGlobal(hit.localPosition());
0243     } else {
0244       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0245     }
0246     LocalPoint lp = refcscchamber()->toLocal(gp);
0247 
0248     // LogDebug
0249     std::stringstream lpss;
0250     lpss << lp;
0251     std::string lps = lpss.str();
0252     std::stringstream gpss;
0253     gpss << gp;
0254     std::string gps = gpss.str();
0255     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit global position " << std::setw(30)
0256                                      << gps << " and local position " << std::setw(30) << lps
0257                                      << " wrt reference csc chamber " << refcscchamber()->id().rawId() << " = "
0258                                      << refcscchamber()->id();
0259 
0260     // Local position of hit w.r.t. chamber
0261     double u = lp.x();
0262     double v = lp.y();
0263     double z = lp.z();
0264 
0265     // Covariance matrix of local errors
0266     SMatrixSym2 IC;  // 2x2, init to 0
0267 
0268     IC(0, 0) = hit.localPositionError().xx();
0269     IC(1, 1) = hit.localPositionError().yy();
0270     //@@ NOT SURE WHICH OFF-DIAGONAL ELEMENT MUST BE DEFINED BUT (1,0) WORKS
0271     //@@ (and SMatrix enforces symmetry)
0272     IC(1, 0) = hit.localPositionError().xy();
0273     // IC(0,1) = IC(1,0);
0274 
0275     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] 2x2 covariance matrix for this Tracking RecHit :: [["
0276                                      << IC(0, 0) << ", " << IC(0, 1) << "][" << IC(1, 0) << "," << IC(1, 1) << "]]";
0277 
0278     // Invert covariance matrix (and trap if it fails!)
0279     bool ok = IC.Invert();
0280     if (!ok) {
0281       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] Failed to invert covariance matrix: \n" << IC;
0282       return;  // MATRIX INVERSION FAILED ... QUIT VOID FUNCTION
0283     }
0284 
0285     // M = (AT E A)
0286     // B = (AT E m)
0287     // for now fill only with sum of blockdiagonal
0288     // elements of (E^-1)_i = IC_i for hit i
0289     M(0, 0) += IC(0, 0);
0290     M(0, 1) += IC(0, 1);
0291     M(0, 2) += IC(0, 0) * z;
0292     M(0, 3) += IC(0, 1) * z;
0293     B(0) += u * IC(0, 0) + v * IC(0, 1);
0294 
0295     M(1, 0) += IC(1, 0);
0296     M(1, 1) += IC(1, 1);
0297     M(1, 2) += IC(1, 0) * z;
0298     M(1, 3) += IC(1, 1) * z;
0299     B(1) += u * IC(1, 0) + v * IC(1, 1);
0300 
0301     M(2, 0) += IC(0, 0) * z;
0302     M(2, 1) += IC(0, 1) * z;
0303     M(2, 2) += IC(0, 0) * z * z;
0304     M(2, 3) += IC(0, 1) * z * z;
0305     B(2) += (u * IC(0, 0) + v * IC(0, 1)) * z;
0306 
0307     M(3, 0) += IC(1, 0) * z;
0308     M(3, 1) += IC(1, 1) * z;
0309     M(3, 2) += IC(1, 0) * z * z;
0310     M(3, 3) += IC(1, 1) * z * z;
0311     B(3) += (u * IC(1, 0) + v * IC(1, 1)) * z;
0312 
0313   }  // End Loop over the TrackingRecHits to make the block matrices to be filled in M and B
0314 
0315   SVector4 p;
0316   bool ok = M.Invert();
0317   if (!ok) {
0318     edm::LogVerbatim("GEMCSCSegment|GEMCSCSegFit") << "[GEMCSCSegFit::fit] Failed to invert matrix: \n" << M;
0319     return;  // MATRIX INVERSION FAILED ... QUIT VOID FUNCTION
0320   } else {
0321     p = M * B;
0322   }
0323 
0324   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] p = " << p(0) << ", " << p(1) << ", " << p(2) << ", "
0325                                    << p(3);
0326 
0327   // fill member variables  (note origin has local z = 0)
0328   //  intercept_
0329   intercept_ = LocalPoint(p(0), p(1), 0.);
0330 
0331   // localdir_ - set so segment points outwards from IP
0332   uslope_ = p(2);
0333   vslope_ = p(3);
0334   setOutFromIP();
0335 
0336   // calculate chi2 of fit
0337   setChi2();
0338 
0339   // flag fit has been done
0340   fitdone_ = true;
0341 }
0342 
0343 void GEMCSCSegFit::setChi2(void) {
0344   double chsq = 0.;
0345   bool gem = false;
0346 
0347   std::vector<const TrackingRecHit*>::const_iterator ih;
0348   for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0349     const TrackingRecHit& hit = (**ih);
0350     GlobalPoint gp;
0351     DetId d = DetId(hit.rawId());
0352     if (d.subdetId() == MuonSubdetId::GEM) {
0353       const GEMEtaPartition* roll = gemetapartition(GEMDetId(d));
0354       gp = roll->toGlobal(hit.localPosition());
0355       gem = true;
0356     } else if (d.subdetId() == MuonSubdetId::CSC) {
0357       const CSCLayer* layer = cscchamber(CSCDetId(d))->layer(CSCDetId(d).layer());
0358       gp = layer->toGlobal(hit.localPosition());
0359       gem = false;
0360     } else {
0361       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0362     }
0363     LocalPoint lp = refcscchamber()->toLocal(gp);
0364 
0365     // LogDebug
0366     std::stringstream lpss;
0367     lpss << lp;
0368     std::string lps = lpss.str();
0369     std::stringstream gpss;
0370     gpss << gp;
0371     std::string gps = gpss.str();
0372     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] - Tracking RecHit in " << (gem ? "GEM" : "CSC")
0373                                      << " global position " << std::setw(30) << gps << " and local position "
0374                                      << std::setw(30) << lps << " wrt reference csc chamber "
0375                                      << refcscchamber()->id().rawId() << " = " << refcscchamber()->id();
0376 
0377     double u = lp.x();
0378     double v = lp.y();
0379     double z = lp.z();
0380 
0381     double du = intercept_.x() + uslope_ * z - u;
0382     double dv = intercept_.y() + vslope_ * z - v;
0383 
0384     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] u, v, z = " << u << ", " << v << ", " << z;
0385 
0386     SMatrixSym2 IC;  // 2x2, init to 0
0387 
0388     IC(0, 0) = hit.localPositionError().xx();
0389     //    IC(0,1) = hit.localPositionError().xy();
0390     IC(1, 0) = hit.localPositionError().xy();
0391     IC(1, 1) = hit.localPositionError().yy();
0392     //    IC(1,0) = IC(0,1);
0393 
0394     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] IC before = \n" << IC;
0395 
0396     // Invert covariance matrix
0397     bool ok = IC.Invert();
0398     if (!ok) {
0399       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Failed to invert covariance matrix: \n" << IC;
0400       return;  // MATRIX INVERSION FAILED ... QUIT VOID FUNCTION
0401     }
0402     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] IC after = \n" << IC;
0403     chsq += du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0404     // LogDebug
0405     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Contribution of this Tracking RecHit to Chi2: "
0406                                         "du^2*D(1,1) + 2*du*dv*D(1,2) + dv^2*D(2,2) = "
0407                                      << du * du << "*" << IC(0, 0) << " + 2.*" << du << "*" << dv << "*" << IC(0, 1)
0408                                      << " + " << dv * dv << "*" << IC(1, 1) << " = "
0409                                      << du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0410   }
0411   // LogDebug
0412   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Total Chi2 = " << chsq;
0413   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Total NDof = " << 2. * hits_.size() - 4;
0414   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Total Chi2/NDof = "
0415                                    << ((hits_.size() > 2) ? (chsq / (2. * hits_.size() - 4)) : (0.0));
0416 
0417   // fill member variables
0418   chi2_ = chsq;
0419   ndof_ = 2. * hits_.size() - 4;
0420 
0421   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] chi2 = " << chi2_ << "/" << ndof_ << " dof";
0422 }
0423 
0424 GEMCSCSegFit::SMatrixSym16 GEMCSCSegFit::weightMatrix() {
0425   bool ok = true;
0426 
0427   SMatrixSym16 matrix = ROOT::Math::SMatrixIdentity();
0428   // for CSC segment ::     max 6 rechits => 12x12,
0429   // for GEM-CSC segment :: max 8 rechits => 16x16
0430   // for all :: init to 1's on diag
0431 
0432   int row = 0;
0433 
0434   for (std::vector<const TrackingRecHit*>::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
0435     const TrackingRecHit& hit = (**it);
0436 
0437     // Note scaleXError allows rescaling the x error if necessary
0438 
0439     matrix(row, row) = scaleXError() * hit.localPositionError().xx();
0440     matrix(row, row + 1) = hit.localPositionError().xy();
0441     ++row;
0442     matrix(row, row - 1) = hit.localPositionError().xy();
0443     matrix(row, row) = hit.localPositionError().yy();
0444     ++row;
0445   }
0446 
0447   ok = matrix.Invert();  // invert in place
0448   if (!ok) {
0449     edm::LogVerbatim("GEMCSCSegment|GEMCSCSegFit") << "[GEMCSCSegFit::weightMatrix] Failed to invert matrix: \n"
0450                                                    << matrix;
0451     SMatrixSym16 emptymatrix = ROOT::Math::SMatrixIdentity();
0452     return emptymatrix;  // return (empty) identity matrix if matrix inversion failed
0453   }
0454   return matrix;
0455 }
0456 
0457 GEMCSCSegFit::SMatrix16by4 GEMCSCSegFit::derivativeMatrix() {
0458   SMatrix16by4 matrix;  // 16x4, init to 0
0459   int row = 0;
0460 
0461   for (std::vector<const TrackingRecHit*>::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
0462     const TrackingRecHit& hit = (**it);
0463     GlobalPoint gp;
0464     DetId d = DetId(hit.rawId());
0465     if (d.subdetId() == MuonSubdetId::GEM) {
0466       const GEMEtaPartition* roll = gemetapartition(GEMDetId(d));
0467       gp = roll->toGlobal(hit.localPosition());
0468     } else if (d.subdetId() == MuonSubdetId::CSC) {
0469       const CSCLayer* layer = cscchamber(CSCDetId(d))->layer(CSCDetId(d).layer());
0470       gp = layer->toGlobal(hit.localPosition());
0471     } else {
0472       edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
0473     }
0474     LocalPoint lp = refcscchamber()->toLocal(gp);
0475     float z = lp.z();
0476 
0477     matrix(row, 0) = 1.;
0478     matrix(row, 2) = z;
0479     ++row;
0480     matrix(row, 1) = 1.;
0481     matrix(row, 3) = z;
0482     ++row;
0483   }
0484   return matrix;
0485 }
0486 
0487 void GEMCSCSegFit::setOutFromIP() {
0488   // Set direction of segment to point from IP outwards
0489   // (Incorrect for particles not coming from IP, of course.)
0490 
0491   double dxdz = uslope_;
0492   double dydz = vslope_;
0493   double dz = 1. / sqrt(1. + dxdz * dxdz + dydz * dydz);
0494   double dx = dz * dxdz;
0495   double dy = dz * dydz;
0496   LocalVector localDir(dx, dy, dz);
0497 
0498   // localDir sometimes needs a sign flip
0499   // Examine its direction and origin in global z: to point outward
0500   // the localDir should always have same sign as global z...
0501 
0502   double globalZpos = (refcscchamber()->toGlobal(intercept_)).z();
0503   double globalZdir = (refcscchamber()->toGlobal(localDir)).z();
0504   double directionSign = globalZpos * globalZdir;
0505   localdir_ = (directionSign * localDir).unit();
0506 }
0507 
0508 AlgebraicSymMatrix GEMCSCSegFit::covarianceMatrix() {
0509   SMatrixSym16 weights = weightMatrix();
0510   SMatrix16by4 A = derivativeMatrix();
0511   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] weights matrix W: \n" << weights;
0512   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] derivatives matrix A: \n" << A;
0513 
0514   // (AT W A)^-1
0515   // e.g. See http://www.phys.ufl.edu/~avery/fitting.html, part I
0516 
0517   bool ok;
0518   SMatrixSym4 result = ROOT::Math::SimilarityT(A, weights);
0519   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] (AT W A): \n" << result;
0520   ok = result.Invert();  // inverts in place
0521   if (!ok) {
0522     edm::LogVerbatim("GEMCSCSegment|GEMCSCSegFit") << "[GEMCSCSegFit::calculateError] Failed to invert matrix: \n"
0523                                                    << result;
0524     AlgebraicSymMatrix emptymatrix(4, 0.);
0525     return emptymatrix;  // return empty matrix if matrix inversion failed
0526   }
0527   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] (AT W A)^-1: \n" << result;
0528 
0529   // reorder components to match TrackingRecHit interface (GEMCSCSegment isa TrackingRecHit)
0530   // i.e. slopes first, then positions
0531   AlgebraicSymMatrix flipped = flipErrors(result);
0532 
0533   return flipped;
0534 }
0535 
0536 AlgebraicSymMatrix GEMCSCSegFit::flipErrors(const SMatrixSym4& a) {
0537   // The GEMCSCSegment needs the error matrix re-arranged to match
0538   // parameters in order (uz, vz, u0, v0)
0539   // where uz, vz = slopes, u0, v0 = intercepts
0540 
0541   edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::flipErrors] input: \n" << a;
0542 
0543   AlgebraicSymMatrix hold(4, 0.);
0544 
0545   for (short j = 0; j != 4; ++j) {
0546     for (short i = 0; i != 4; ++i) {
0547       hold(i + 1, j + 1) = a(i, j);  // SMatrix counts from 0, AlgebraicMatrix from 1
0548     }
0549   }
0550 
0551   //  LogTrace("GEMCSCSegFit") << "[GEMCSCSegFit::flipErrors] after copy:";
0552   //  LogTrace("GEMCSCSegFit") << "(" << hold(1,1) << "  " << hold(1,2) << "  " << hold(1,3) << "  " << hold(1,4);
0553   //  LogTrace("GEMCSCSegFit") << " " << hold(2,1) << "  " << hold(2,2) << "  " << hold(2,3) << "  " << hold(2,4);
0554   //  LogTrace("GEMCSCSegFit") << " " << hold(3,1) << "  " << hold(3,2) << "  " << hold(3,3) << "  " << hold(3,4);
0555   //  LogTrace("GEMCSCSegFit") << " " << hold(4,1) << "  " << hold(4,2) << "  " << hold(4,3) << "  " << hold(4,4) << ")";
0556 
0557   // errors on slopes into upper left
0558   hold(1, 1) = a(2, 2);
0559   hold(1, 2) = a(2, 3);
0560   hold(2, 1) = a(3, 2);
0561   hold(2, 2) = a(3, 3);
0562 
0563   // errors on positions into lower right
0564   hold(3, 3) = a(0, 0);
0565   hold(3, 4) = a(0, 1);
0566   hold(4, 3) = a(1, 0);
0567   hold(4, 4) = a(1, 1);
0568 
0569   // must also interchange off-diagonal elements of off-diagonal 2x2 submatrices
0570   hold(4, 1) = a(1, 2);
0571   hold(3, 2) = a(0, 3);
0572   hold(2, 3) = a(3, 0);  // = a(0,3)
0573   hold(1, 4) = a(2, 1);  // = a(1,2)
0574 
0575   //  LogTrace("GEMCSCSegFit") << "[GEMCSCSegFit::flipErrors] after flip:";
0576   //  LogTrace("GEMCSCSegFit") << "(" << hold(1,1) << "  " << hold(1,2) << "  " << hold(1,3) << "  " << hold(1,4);
0577   //  LogTrace("GEMCSCSegFit") << " " << hold(2,1) << "  " << hold(2,2) << "  " << hold(2,3) << "  " << hold(2,4);
0578   //  LogTrace("GEMCSCSegFit") << " " << hold(3,1) << "  " << hold(3,2) << "  " << hold(3,3) << "  " << hold(3,4);
0579   //  LogTrace("GEMCSCSegFit") << " " << hold(4,1) << "  " << hold(4,2) << "  " << hold(4,3) << "  " << hold(4,4) << ")";
0580 
0581   return hold;
0582 }
0583 
0584 float GEMCSCSegFit::xfit(float z) const {
0585   //@@ ADD THIS TO EACH ACCESSOR OF FIT RESULTS?
0586   //  if ( !fitdone() ) fit();
0587   return intercept_.x() + uslope_ * z;
0588 }
0589 
0590 float GEMCSCSegFit::yfit(float z) const { return intercept_.y() + vslope_ * z; }
0591 
0592 float GEMCSCSegFit::xdev(float x, float z) const { return intercept_.x() + uslope_ * z - x; }
0593 
0594 float GEMCSCSegFit::ydev(float y, float z) const { return intercept_.y() + vslope_ * z - y; }
0595 
0596 float GEMCSCSegFit::Rdev(float x, float y, float z) const {
0597   return sqrt(xdev(x, z) * xdev(x, z) + ydev(y, z) * ydev(y, z));
0598 }
[ 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!