CMS 3D CMS Logo

GEMCSCSegFit.cc
Go to the documentation of this file.
1 // ------------------------- //
2 // --> GEMCSCSegFit.cc
3 // Created: 21.04.2015
4 // --> Based on CSCSegFit.cc
5 // with last mod: 03.02.2015
6 // as found in 750pre2 rel
7 // ------------------------- //
8 
12 
15 
18 
19 void GEMCSCSegFit::fit(void) {
20  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] - start the fitting fun :: nhits = " << nhits();
21  if (fitdone())
22  return; // don't redo fit unnecessarily
23  short n = nhits();
24  switch (n) {
25  case 1:
26  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] - cannot fit just 1 hit!!";
27  break;
28  case 2:
29  fit2();
30  break;
31  case 3:
32  case 4:
33  case 5:
34  case 6:
35  case 7:
36  case 8:
37  fitlsq();
38  break;
39  default:
40  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] - cannot fit more than 8 hits!!";
41  }
42 }
43 
44 void GEMCSCSegFit::fit2(void) {
45  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit2] - start fit2()";
46  // Just join the two points
47  // Equation of straight line between (x1, y1) and (x2, y2) in xy-plane is
48  // y = mx + c
49  // with m = (y2-y1)/(x2-x1)
50  // and c = (y1*x2-x2*y1)/(x2-x1)
51 
52  // 1) Check whether hits are on the same layer
53  // -------------------------------------------
54  std::vector<const TrackingRecHit*>::const_iterator ih = hits_.begin();
55  // layer numbering: GEM: (1,2) CSC (3,4,5,6,7,8)
56  int il1 = 0, il2 = 0;
57  DetId d1 = DetId((*ih)->rawId());
58  // check whether first hit is GEM or CSC
59  if (d1.subdetId() == MuonSubdetId::GEM) {
60  il1 = GEMDetId(d1).layer();
61  } else if (d1.subdetId() == MuonSubdetId::CSC) {
62  il1 = CSCDetId(d1).layer() + 2;
63  } else {
64  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
65  }
66  const TrackingRecHit& h1 = (**ih);
67  ++ih;
68  DetId d2 = DetId((*ih)->rawId());
69  // check whether second hit is GEM or CSC
70  if (d2.subdetId() == MuonSubdetId::GEM) {
71  il2 = GEMDetId(d2).layer();
72  } else if (d2.subdetId() == MuonSubdetId::CSC) {
73  il2 = GEMDetId(d2).layer() + 2;
74  } else {
75  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
76  }
77  const TrackingRecHit& h2 = (**ih);
78  // Skip if on same layer, but should be impossible :)
79  if (il1 == il2) {
80  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - 2 hits on same layer!!";
81  return;
82  }
83 
84  // 2) Global Positions of hit 1 and 2 and
85  // Local Positions of hit 1 and 2 w.r.t. reference CSC Chamber Frame
86  // ---------------------------------------------------------------------
87  GlobalPoint h1glopos, h2glopos;
88  // global position hit 1
89  if (d1.subdetId() == MuonSubdetId::GEM) {
90  const GEMEtaPartition* roll1 = gemetapartition(GEMDetId(d1));
91  h1glopos = roll1->toGlobal(h1.localPosition());
92  } else if (d1.subdetId() == MuonSubdetId::CSC) {
93  const CSCLayer* layer1 = cscchamber(CSCDetId(d1))->layer(CSCDetId(d1).layer());
94  h1glopos = layer1->toGlobal(h1.localPosition());
95  } else {
96  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
97  }
98  // global position hit 2
99  if (d2.subdetId() == MuonSubdetId::GEM) {
100  const GEMEtaPartition* roll2 = gemetapartition(GEMDetId(d2));
101  h2glopos = roll2->toGlobal(h2.localPosition());
102  } else if (d2.subdetId() == MuonSubdetId::CSC) {
103  const CSCLayer* layer2 = cscchamber(CSCDetId(d2))->layer(CSCDetId(d2).layer());
104  h2glopos = layer2->toGlobal(h2.localPosition());
105  } else {
106  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
107  }
108  // local positions hit 1 and 2 w.r.t. ref CSC Chamber Frame
109  // We want hit wrt chamber (and local z will be != 0)
110  LocalPoint h1pos = refcscchamber()->toLocal(h1glopos);
111  LocalPoint h2pos = refcscchamber()->toLocal(h2glopos);
112 
113  // 3) Now make straight line between the two points in local coords
114  // ----------------------------------------------------------------
115  float dz = h2pos.z() - h1pos.z();
116  if (dz != 0) {
117  uslope_ = (h2pos.x() - h1pos.x()) / dz;
118  vslope_ = (h2pos.y() - h1pos.y()) / dz;
119  }
120  float uintercept = (h1pos.x() * h2pos.z() - h2pos.x() * h1pos.z()) / dz;
121  float vintercept = (h1pos.y() * h2pos.z() - h2pos.y() * h1pos.z()) / dz;
122  intercept_ = LocalPoint(uintercept, vintercept, 0.);
123 
124  setOutFromIP();
125 
126  //@@ NOT SURE WHAT IS SENSIBLE FOR THESE...
127  chi2_ = 0.;
128  ndof_ = 0;
129 
130  fitdone_ = true;
131 }
132 
134  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - start fitlsq";
135  // Linear least-squares fit to up to 6 CSC rechits, one per layer in a CSC
136  // and up to 2 GEM rechits in a GEM superchamber
137  // (we can later later on go up to 4 hits in case of an overlapping chamber)
138  // (... and if there is an overlap in the GEM chambers, than maybe we also
139  // have an overlap in the CSC chambers... maybe we can also benefit there)
140 
141  // Comments below taken from CSCSegFit algorithm.
142  // Comments adapted from original CSCSegAlgoSK algorithm.
143 
144  // Fit to the local x, y rechit coordinates in z projection
145  // The CSC & GEM strip measurement controls the precision of x
146  // The CSC wire measurement & GEM eta-partition controls the precision of y.
147  // Typical precision CSC: u (strip, sigma~200um), v (wire, sigma~1cm)
148  // Typical precision GEM: u (strip, sigma~250um), v (eta-part, sigma~10-20cm)
149 
150  // Set up the normal equations for the least-squares fit as a matrix equation
151 
152  // We have a vector of measurements m, which is a 2n x 1 dim matrix
153  // The transpose mT is (u1, v1, u2, v2, ..., un, vn) where
154  // ui is the strip-associated measurement and
155  // vi is the wire-associated measurement
156  // for a given rechit i.
157 
158  // The fit is to
159  // u = u0 + uz * z
160  // v = v0 + vz * z
161  // where u0, uz, v0, vz are the parameters to be obtained from the fit.
162 
163  // These are contained in a vector p which is a 4x1 dim matrix, and
164  // its transpose pT is (u0, v0, uz, vz). Note the ordering!
165 
166  // The covariance matrix for each pair of measurements is 2 x 2 and
167  // the inverse of this is the error matrix E.
168  // The error matrix for the whole set of n measurements is a diagonal
169  // matrix with diagonal elements the individual 2 x 2 error matrices
170  // (because the inverse of a diagonal matrix is a diagonal matrix
171  // with each element the inverse of the original.)
172 
173  // In function 'weightMatrix()', the variable 'matrix' is filled with this
174  // block-diagonal overall covariance matrix. Then 'matrix' is inverted to the
175  // block-diagonal error matrix, and returned.
176 
177  // Define the matrix A as
178  // 1 0 z1 0
179  // 0 1 0 z1
180  // 1 0 z2 0
181  // 0 1 0 z2
182  // .. .. .. ..
183  // 1 0 zn 0
184  // 0 1 0 zn
185 
186  // This matrix A is set up and returned by function 'derivativeMatrix()'.
187 
188  // Then the normal equations are described by the matrix equation
189  //
190  // (AT E A)p = (AT E)m
191  //
192  // where AT is the transpose of A.
193 
194  // Call the combined matrix on the LHS, M, and that on the RHS, B:
195  // M p = B
196 
197  // We solve this for the parameter vector, p.
198  // The elements of M and B then involve sums over the hits
199 
200  // The covariance matrix of the parameters is obtained by
201  // (AT E A)^-1 calculated in 'covarianceMatrix()'.
202 
203  // NOTE
204  // We need local position of a RecHit w.r.t. the CHAMBER
205  // and the RecHit itself only knows its local position w.r.t.
206  // the LAYER, so we must explicitly transform global position.
207 
208  SMatrix4 M; // 4x4, init to 0
209  SVector4 B; // 4x1, init to 0;
210 
211  std::vector<const TrackingRecHit*>::const_iterator ih = hits_.begin();
212 
213  // LogDebug :: Loop over the TrackingRecHits and print the GEM and CSC Hits
214  for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
215  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - looping over TrackingRecHits";
216  const TrackingRecHit& hit = (**ih);
217  DetId d = DetId(hit.rawId());
218  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit in detid (" << d.rawId() << ")";
219  if (d.subdetId() == MuonSubdetId::GEM) {
220  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit is a GEM Hit in detid ("
221  << d.rawId() << ")";
222  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - GEM DetId (" << GEMDetId(d.rawId()) << ")";
223  } else if (d.subdetId() == MuonSubdetId::CSC) {
224  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit is a CSC Hit in detid ("
225  << d.rawId() << ")";
226  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - CSC DetId (" << CSCDetId(d.rawId()) << ")";
227  } else {
228  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
229  }
230  }
231 
232  // Loop over the TrackingRecHits and make small (2x2) matrices used to fill the blockdiagonal covariance matrix E^-1
233  for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
234  const TrackingRecHit& hit = (**ih);
235  GlobalPoint gp;
236  DetId d = DetId(hit.rawId());
237  if (d.subdetId() == MuonSubdetId::GEM) {
239  gp = roll->toGlobal(hit.localPosition());
240  } else if (d.subdetId() == MuonSubdetId::CSC) {
242  gp = layer->toGlobal(hit.localPosition());
243  } else {
244  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
245  }
247 
248  // LogDebug
249  std::stringstream lpss;
250  lpss << lp;
251  std::string lps = lpss.str();
252  std::stringstream gpss;
253  gpss << gp;
254  std::string gps = gpss.str();
255  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fitlsq] - Tracking RecHit global position " << std::setw(30)
256  << gps << " and local position " << std::setw(30) << lps
257  << " wrt reference csc chamber " << refcscchamber()->id().rawId() << " = "
258  << refcscchamber()->id();
259 
260  // Local position of hit w.r.t. chamber
261  double u = lp.x();
262  double v = lp.y();
263  double z = lp.z();
264 
265  // Covariance matrix of local errors
266  SMatrixSym2 IC; // 2x2, init to 0
267 
268  IC(0, 0) = hit.localPositionError().xx();
269  IC(1, 1) = hit.localPositionError().yy();
270  //@@ NOT SURE WHICH OFF-DIAGONAL ELEMENT MUST BE DEFINED BUT (1,0) WORKS
271  //@@ (and SMatrix enforces symmetry)
272  IC(1, 0) = hit.localPositionError().xy();
273  // IC(0,1) = IC(1,0);
274 
275  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] 2x2 covariance matrix for this Tracking RecHit :: [["
276  << IC(0, 0) << ", " << IC(0, 1) << "][" << IC(1, 0) << "," << IC(1, 1) << "]]";
277 
278  // Invert covariance matrix (and trap if it fails!)
279  bool ok = IC.Invert();
280  if (!ok) {
281  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] Failed to invert covariance matrix: \n" << IC;
282  return; // MATRIX INVERSION FAILED ... QUIT VOID FUNCTION
283  }
284 
285  // M = (AT E A)
286  // B = (AT E m)
287  // for now fill only with sum of blockdiagonal
288  // elements of (E^-1)_i = IC_i for hit i
289  M(0, 0) += IC(0, 0);
290  M(0, 1) += IC(0, 1);
291  M(0, 2) += IC(0, 0) * z;
292  M(0, 3) += IC(0, 1) * z;
293  B(0) += u * IC(0, 0) + v * IC(0, 1);
294 
295  M(1, 0) += IC(1, 0);
296  M(1, 1) += IC(1, 1);
297  M(1, 2) += IC(1, 0) * z;
298  M(1, 3) += IC(1, 1) * z;
299  B(1) += u * IC(1, 0) + v * IC(1, 1);
300 
301  M(2, 0) += IC(0, 0) * z;
302  M(2, 1) += IC(0, 1) * z;
303  M(2, 2) += IC(0, 0) * z * z;
304  M(2, 3) += IC(0, 1) * z * z;
305  B(2) += (u * IC(0, 0) + v * IC(0, 1)) * z;
306 
307  M(3, 0) += IC(1, 0) * z;
308  M(3, 1) += IC(1, 1) * z;
309  M(3, 2) += IC(1, 0) * z * z;
310  M(3, 3) += IC(1, 1) * z * z;
311  B(3) += (u * IC(1, 0) + v * IC(1, 1)) * z;
312 
313  } // End Loop over the TrackingRecHits to make the block matrices to be filled in M and B
314 
315  SVector4 p;
316  bool ok = M.Invert();
317  if (!ok) {
318  edm::LogVerbatim("GEMCSCSegment|GEMCSCSegFit") << "[GEMCSCSegFit::fit] Failed to invert matrix: \n" << M;
319  return; // MATRIX INVERSION FAILED ... QUIT VOID FUNCTION
320  } else {
321  p = M * B;
322  }
323 
324  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::fit] p = " << p(0) << ", " << p(1) << ", " << p(2) << ", "
325  << p(3);
326 
327  // fill member variables (note origin has local z = 0)
328  // intercept_
329  intercept_ = LocalPoint(p(0), p(1), 0.);
330 
331  // localdir_ - set so segment points outwards from IP
332  uslope_ = p(2);
333  vslope_ = p(3);
334  setOutFromIP();
335 
336  // calculate chi2 of fit
337  setChi2();
338 
339  // flag fit has been done
340  fitdone_ = true;
341 }
342 
343 void GEMCSCSegFit::setChi2(void) {
344  double chsq = 0.;
345  bool gem = false;
346 
347  std::vector<const TrackingRecHit*>::const_iterator ih;
348  for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
349  const TrackingRecHit& hit = (**ih);
350  GlobalPoint gp;
351  DetId d = DetId(hit.rawId());
352  if (d.subdetId() == MuonSubdetId::GEM) {
354  gp = roll->toGlobal(hit.localPosition());
355  gem = true;
356  } else if (d.subdetId() == MuonSubdetId::CSC) {
358  gp = layer->toGlobal(hit.localPosition());
359  gem = false;
360  } else {
361  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
362  }
364 
365  // LogDebug
366  std::stringstream lpss;
367  lpss << lp;
368  std::string lps = lpss.str();
369  std::stringstream gpss;
370  gpss << gp;
371  std::string gps = gpss.str();
372  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] - Tracking RecHit in " << (gem ? "GEM" : "CSC")
373  << " global position " << std::setw(30) << gps << " and local position "
374  << std::setw(30) << lps << " wrt reference csc chamber "
375  << refcscchamber()->id().rawId() << " = " << refcscchamber()->id();
376 
377  double u = lp.x();
378  double v = lp.y();
379  double z = lp.z();
380 
381  double du = intercept_.x() + uslope_ * z - u;
382  double dv = intercept_.y() + vslope_ * z - v;
383 
384  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] u, v, z = " << u << ", " << v << ", " << z;
385 
386  SMatrixSym2 IC; // 2x2, init to 0
387 
388  IC(0, 0) = hit.localPositionError().xx();
389  // IC(0,1) = hit.localPositionError().xy();
390  IC(1, 0) = hit.localPositionError().xy();
391  IC(1, 1) = hit.localPositionError().yy();
392  // IC(1,0) = IC(0,1);
393 
394  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] IC before = \n" << IC;
395 
396  // Invert covariance matrix
397  bool ok = IC.Invert();
398  if (!ok) {
399  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Failed to invert covariance matrix: \n" << IC;
400  return; // MATRIX INVERSION FAILED ... QUIT VOID FUNCTION
401  }
402  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] IC after = \n" << IC;
403  chsq += du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
404  // LogDebug
405  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Contribution of this Tracking RecHit to Chi2: "
406  "du^2*D(1,1) + 2*du*dv*D(1,2) + dv^2*D(2,2) = "
407  << du * du << "*" << IC(0, 0) << " + 2.*" << du << "*" << dv << "*" << IC(0, 1)
408  << " + " << dv * dv << "*" << IC(1, 1) << " = "
409  << du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
410  }
411  // LogDebug
412  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Total Chi2 = " << chsq;
413  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Total NDof = " << 2. * hits_.size() - 4;
414  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] Total Chi2/NDof = "
415  << ((hits_.size() > 2) ? (chsq / (2. * hits_.size() - 4)) : (0.0));
416 
417  // fill member variables
418  chi2_ = chsq;
419  ndof_ = 2. * hits_.size() - 4;
420 
421  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::setChi2] chi2 = " << chi2_ << "/" << ndof_ << " dof";
422 }
423 
425  bool ok = true;
426 
427  SMatrixSym16 matrix = ROOT::Math::SMatrixIdentity();
428  // for CSC segment :: max 6 rechits => 12x12,
429  // for GEM-CSC segment :: max 8 rechits => 16x16
430  // for all :: init to 1's on diag
431 
432  int row = 0;
433 
434  for (std::vector<const TrackingRecHit*>::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
435  const TrackingRecHit& hit = (**it);
436 
437  // Note scaleXError allows rescaling the x error if necessary
438 
439  matrix(row, row) = scaleXError() * hit.localPositionError().xx();
440  matrix(row, row + 1) = hit.localPositionError().xy();
441  ++row;
442  matrix(row, row - 1) = hit.localPositionError().xy();
443  matrix(row, row) = hit.localPositionError().yy();
444  ++row;
445  }
446 
447  ok = matrix.Invert(); // invert in place
448  if (!ok) {
449  edm::LogVerbatim("GEMCSCSegment|GEMCSCSegFit") << "[GEMCSCSegFit::weightMatrix] Failed to invert matrix: \n"
450  << matrix;
451  SMatrixSym16 emptymatrix = ROOT::Math::SMatrixIdentity();
452  return emptymatrix; // return (empty) identity matrix if matrix inversion failed
453  }
454  return matrix;
455 }
456 
458  SMatrix16by4 matrix; // 16x4, init to 0
459  int row = 0;
460 
461  for (std::vector<const TrackingRecHit*>::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
462  const TrackingRecHit& hit = (**it);
463  GlobalPoint gp;
464  DetId d = DetId(hit.rawId());
465  if (d.subdetId() == MuonSubdetId::GEM) {
467  gp = roll->toGlobal(hit.localPosition());
468  } else if (d.subdetId() == MuonSubdetId::CSC) {
470  gp = layer->toGlobal(hit.localPosition());
471  } else {
472  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit:fit2] - TrackingRecHit is not in GEM or CSC subdetector";
473  }
475  float z = lp.z();
476 
477  matrix(row, 0) = 1.;
478  matrix(row, 2) = z;
479  ++row;
480  matrix(row, 1) = 1.;
481  matrix(row, 3) = z;
482  ++row;
483  }
484  return matrix;
485 }
486 
488  // Set direction of segment to point from IP outwards
489  // (Incorrect for particles not coming from IP, of course.)
490 
491  double dxdz = uslope_;
492  double dydz = vslope_;
493  double dz = 1. / sqrt(1. + dxdz * dxdz + dydz * dydz);
494  double dx = dz * dxdz;
495  double dy = dz * dydz;
496  LocalVector localDir(dx, dy, dz);
497 
498  // localDir sometimes needs a sign flip
499  // Examine its direction and origin in global z: to point outward
500  // the localDir should always have same sign as global z...
501 
502  double globalZpos = (refcscchamber()->toGlobal(intercept_)).z();
503  double globalZdir = (refcscchamber()->toGlobal(localDir)).z();
504  double directionSign = globalZpos * globalZdir;
505  localdir_ = (directionSign * localDir).unit();
506 }
507 
511  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] weights matrix W: \n" << weights;
512  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] derivatives matrix A: \n" << A;
513 
514  // (AT W A)^-1
515  // e.g. See http://www.phys.ufl.edu/~avery/fitting.html, part I
516 
517  bool ok;
518  SMatrixSym4 result = ROOT::Math::SimilarityT(A, weights);
519  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] (AT W A): \n" << result;
520  ok = result.Invert(); // inverts in place
521  if (!ok) {
522  edm::LogVerbatim("GEMCSCSegment|GEMCSCSegFit") << "[GEMCSCSegFit::calculateError] Failed to invert matrix: \n"
523  << result;
524  AlgebraicSymMatrix emptymatrix(4, 0.);
525  return emptymatrix; // return empty matrix if matrix inversion failed
526  }
527  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::covarianceMatrix] (AT W A)^-1: \n" << result;
528 
529  // reorder components to match TrackingRecHit interface (GEMCSCSegment isa TrackingRecHit)
530  // i.e. slopes first, then positions
532 
533  return flipped;
534 }
535 
536 AlgebraicSymMatrix GEMCSCSegFit::flipErrors(const SMatrixSym4& a) {
537  // The GEMCSCSegment needs the error matrix re-arranged to match
538  // parameters in order (uz, vz, u0, v0)
539  // where uz, vz = slopes, u0, v0 = intercepts
540 
541  edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::flipErrors] input: \n" << a;
542 
543  AlgebraicSymMatrix hold(4, 0.);
544 
545  for (short j = 0; j != 4; ++j) {
546  for (short i = 0; i != 4; ++i) {
547  hold(i + 1, j + 1) = a(i, j); // SMatrix counts from 0, AlgebraicMatrix from 1
548  }
549  }
550 
551  // LogTrace("GEMCSCSegFit") << "[GEMCSCSegFit::flipErrors] after copy:";
552  // LogTrace("GEMCSCSegFit") << "(" << hold(1,1) << " " << hold(1,2) << " " << hold(1,3) << " " << hold(1,4);
553  // LogTrace("GEMCSCSegFit") << " " << hold(2,1) << " " << hold(2,2) << " " << hold(2,3) << " " << hold(2,4);
554  // LogTrace("GEMCSCSegFit") << " " << hold(3,1) << " " << hold(3,2) << " " << hold(3,3) << " " << hold(3,4);
555  // LogTrace("GEMCSCSegFit") << " " << hold(4,1) << " " << hold(4,2) << " " << hold(4,3) << " " << hold(4,4) << ")";
556 
557  // errors on slopes into upper left
558  hold(1, 1) = a(2, 2);
559  hold(1, 2) = a(2, 3);
560  hold(2, 1) = a(3, 2);
561  hold(2, 2) = a(3, 3);
562 
563  // errors on positions into lower right
564  hold(3, 3) = a(0, 0);
565  hold(3, 4) = a(0, 1);
566  hold(4, 3) = a(1, 0);
567  hold(4, 4) = a(1, 1);
568 
569  // must also interchange off-diagonal elements of off-diagonal 2x2 submatrices
570  hold(4, 1) = a(1, 2);
571  hold(3, 2) = a(0, 3);
572  hold(2, 3) = a(3, 0); // = a(0,3)
573  hold(1, 4) = a(2, 1); // = a(1,2)
574 
575  // LogTrace("GEMCSCSegFit") << "[GEMCSCSegFit::flipErrors] after flip:";
576  // LogTrace("GEMCSCSegFit") << "(" << hold(1,1) << " " << hold(1,2) << " " << hold(1,3) << " " << hold(1,4);
577  // LogTrace("GEMCSCSegFit") << " " << hold(2,1) << " " << hold(2,2) << " " << hold(2,3) << " " << hold(2,4);
578  // LogTrace("GEMCSCSegFit") << " " << hold(3,1) << " " << hold(3,2) << " " << hold(3,3) << " " << hold(3,4);
579  // LogTrace("GEMCSCSegFit") << " " << hold(4,1) << " " << hold(4,2) << " " << hold(4,3) << " " << hold(4,4) << ")";
580 
581  return hold;
582 }
583 
584 float GEMCSCSegFit::xfit(float z) const {
585  //@@ ADD THIS TO EACH ACCESSOR OF FIT RESULTS?
586  // if ( !fitdone() ) fit();
587  return intercept_.x() + uslope_ * z;
588 }
589 
590 float GEMCSCSegFit::yfit(float z) const { return intercept_.y() + vslope_ * z; }
591 
592 float GEMCSCSegFit::xdev(float x, float z) const { return intercept_.x() + uslope_ * z - x; }
593 
594 float GEMCSCSegFit::ydev(float y, float z) const { return intercept_.y() + vslope_ * z - y; }
595 
596 float GEMCSCSegFit::Rdev(float x, float y, float z) const {
597  return sqrt(xdev(x, z) * xdev(x, z) + ydev(y, z) * ydev(y, z));
598 }
void setChi2(void)
Log< level::Info, true > LogVerbatim
float dydz
const CSCLayer * layer(CSCDetId id) const
Return the layer corresponding to the given id.
Definition: CSCChamber.cc:30
SMatrixSym16 weightMatrix(void)
std::vector< const TrackingRecHit * > hits_
Definition: GEMCSCSegFit.h:196
Point3DBase< Scalar, LocalTag > LocalPoint
Definition: Definitions.h:30
static constexpr int GEM
Definition: MuonSubdetId.h:14
Definition: APVGainStruct.h:7
float xdev(float x, float z) const
float dxdz
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:58
size_t nhits(void) const
Definition: GEMCSCSegFit.h:133
T z() const
Definition: PV3DBase.h:61
void fitlsq(void)
float ydev(float y, float z) const
SMatrix16by4 derivativeMatrix(void)
void fit(void)
Definition: GEMCSCSegFit.cc:19
int layer() const
Definition: CSCDetId.h:56
ROOT::Math::SVector< double, 4 > SVector4
Definition: GEMCSCSegFit.h:65
constexpr int layer() const
Definition: GEMDetId.h:190
bool fitdone() const
Definition: GEMCSCSegFit.h:171
float yfit(float z) const
const GEMEtaPartition * gemetapartition(uint32_t id) const
Definition: GEMCSCSegFit.h:155
AlgebraicSymMatrix covarianceMatrix(void)
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:23
LocalVector localdir_
Definition: GEMCSCSegFit.h:200
float Rdev(float x, float y, float z) const
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
const CSCChamber * refcscchamber() const
Definition: GEMCSCSegFit.h:163
Basic3DVector unit() const
d
Definition: ztail.py:151
const CSCChamber * cscchamber(uint32_t id) const
Definition: GEMCSCSegFit.h:139
ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2 > > SMatrixSym2
Definition: GEMCSCSegFit.h:62
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
Definition: DetId.h:17
AlgebraicSymMatrix flipErrors(const SMatrixSym4 &)
ROOT::Math::SMatrix< double, 16, 4 > SMatrix16by4
Definition: GEMCSCSegFit.h:55
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
ROOT::Math::SMatrix< double, 4 > SMatrix4
Definition: GEMCSCSegFit.h:58
LocalPoint intercept_
Definition: GEMCSCSegFit.h:199
void setOutFromIP(void)
double a
Definition: hdecay.h:121
CLHEP::HepSymMatrix AlgebraicSymMatrix
float x
ROOT::Math::SMatrix< double, 16, 16, ROOT::Math::MatRepSym< double, 16 > > SMatrixSym16
Definition: GEMCSCSegFit.h:52
Definition: APVGainStruct.h:7
float xfit(float z) const
ROOT::Math::SMatrix< double, 4, 4, ROOT::Math::MatRepSym< double, 4 > > SMatrixSym4
Definition: GEMCSCSegFit.h:59
static constexpr float d1
double scaleXError(void) const
Definition: GEMCSCSegFit.h:132
static constexpr int CSC
Definition: MuonSubdetId.h:12
virtual LocalPoint localPosition() const =0
void fit2(void)
Definition: GEMCSCSegFit.cc:44
CSCDetId id() const
Get the (concrete) DetId.
Definition: CSCChamber.h:34