CMS 3D CMS Logo

MuonResidualsFromTrack.cc
Go to the documentation of this file.
1 /*
2  * $Id: $
3  */
4 
8 
15 
18 
19 #include "TDecompChol.h"
20 #include <cmath>
21 
22 
26  edm::ESHandle<DetIdAssociator> muonDetIdAssociator_,
28  const Trajectory *traj,
29  const reco::Track* recoTrack,
31  double maxResidual)
32 : m_recoTrack(recoTrack)
33 {
34  bool m_debug = false;
35 
36  if (m_debug) {
37  std::cout << "BEGIN MuonResidualsFromTrack" << std::endl;
38  const std::string metname = " *** MuonResidualsFromTrack *** ";
39  LogTrace(metname) << "Tracking Component changed!";
40  }
41 
42  clear();
43 
44  edm::ESHandle<TransientTrackingRecHitBuilder> theTrackerRecHitBuilder;
45  iSetup.get<TransientRecHitRecord>().get("WithTrackAngle",theTrackerRecHitBuilder);
46  reco::TransientTrack track( *m_recoTrack, &*magneticField, globalGeometry );
48  int iT = 0, iM = 0;
50  if((*hit)->isValid()) {
51  DetId hitId = (*hit)->geographicalId();
52  if ( hitId.det() == DetId::Tracker ) {
53  iT++;
54  if (m_debug) std::cout << "Tracker Hit " << iT << " is found. Add to refit. Dimension: " << (*hit)->dimension() << std::endl;
55 
56  recHitsForRefit.push_back( theTrackerRecHitBuilder->build(&**hit) );
57  } else if ( hitId.det() == DetId::Muon ){
58  // if ( (*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit ) {
59  // LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged";
60  // continue;
61  // }
62  iM++;
63  if (m_debug) std::cout << "Muon Hit " << iM << " is found. We do not add muon hits to refit. Dimension: " << (*hit)->dimension() << std::endl;
64  if ( hitId.subdetId() == MuonSubdetId::DT ) {
65  const DTChamberId chamberId(hitId.rawId());
66  if (m_debug) std::cout << "Muon Hit in DT wheel " << chamberId.wheel() << " station " << chamberId.station() << " sector " << chamberId.sector() << "." << std::endl;
67  } else if ( hitId.subdetId() == MuonSubdetId::CSC ) {
68  const CSCDetId cscDetId(hitId.rawId());
69  if (m_debug) std::cout << "Muon hit in CSC endcap " << cscDetId.endcap() << " station " << cscDetId.station() << " ring " << cscDetId.ring() << " chamber " << cscDetId.chamber() << "." << std::endl;
70  } else if ( hitId.subdetId() == MuonSubdetId::RPC ) {
71  if (m_debug) std::cout << "Muon Hit in RPC" << std::endl;
72  } else {
73  if (m_debug) std::cout << "Warning! Muon Hit not in DT or CSC or RPC" << std::endl;
74  }
75  // recHitsForRefit.push_back(theMuonRecHitBuilder->build(&**hit));
76  }
77  }
78  }
79 
80  // TrackTransformer trackTransformer();
81  // std::vector<Trajectory> vTrackerTrajectory = trackTransformer.transform(track, recHitsForReFit);
82  // std::cout << "Tracker trajectories size " << vTrackerTrajectory.size() << std::endl;
83 
84 
85  TrajectoryStateOnSurface lastTrackerTsos;
86  double lastTrackerTsosGlobalPositionR = 0.0;
87 
88  std::vector<TrajectoryMeasurement> vTrajMeasurement = traj->measurements();
89  if (m_debug) std::cout << " Size of vector of TrajectoryMeasurements: " << vTrajMeasurement.size() << std::endl;
90  int nTrajMeasurement = 0;
91  for ( std::vector<TrajectoryMeasurement>::const_iterator iTrajMeasurement = vTrajMeasurement.begin();
92  iTrajMeasurement != vTrajMeasurement.end();
93  ++iTrajMeasurement ) {
94  nTrajMeasurement++;
95  if (m_debug) std::cout << " TrajectoryMeasurement #" << nTrajMeasurement << std::endl;
96 
97  TrajectoryMeasurement trajMeasurement = *iTrajMeasurement;
98 
99  TrajectoryStateOnSurface tsos = m_tsoscomb(trajMeasurement.forwardPredictedState(), trajMeasurement.backwardPredictedState());
100  const TrajectoryStateOnSurface& tsosF = trajMeasurement.forwardPredictedState();
101  const TrajectoryStateOnSurface& tsosB = trajMeasurement.backwardPredictedState();
102  const TrajectoryStateOnSurface& tsosU = trajMeasurement.updatedState();
103  if (m_debug) std::cout << " TrajectoryMeasurement TSOS validity: " << tsos.isValid() << std::endl;
104  if ( tsos.isValid() ) {
105  double tsosGlobalPositionR = sqrt( tsos.globalPosition().x()*tsos.globalPosition().x() + tsos.globalPosition().y()*tsos.globalPosition().y() );
106  if (m_debug) {
107  std::cout << " TrajectoryMeasurement TSOS localPosition"
108  << " x: " << tsos.localPosition().x()
109  << " y: " << tsos.localPosition().y()
110  << " z: " << tsos.localPosition().z()
111  << std::endl;
112  std::cout << " TrajectoryMeasurement TSOS globalPosition"
113  << " x: " << tsos.globalPosition().x()
114  << " y: " << tsos.globalPosition().y()
115  << " R: " << tsosGlobalPositionR
116  << " z: " << tsos.globalPosition().z()
117  << std::endl;
118  }
119  if ( tsosGlobalPositionR > lastTrackerTsosGlobalPositionR ) {
120  lastTrackerTsos = tsos;
121  lastTrackerTsosGlobalPositionR = tsosGlobalPositionR;
122  }
123  }
124 
125  const TransientTrackingRecHit *trajMeasurementHit = &(*trajMeasurement.recHit());
126  if (m_debug) std::cout << " TrajectoryMeasurement hit validity: " << trajMeasurementHit->isValid() << std::endl;
127  if ( trajMeasurementHit->isValid() ) {
128  DetId trajMeasurementHitId = trajMeasurementHit->geographicalId();
129  int trajMeasurementHitDim = trajMeasurementHit->dimension();
130  if ( trajMeasurementHitId.det() == DetId::Tracker ) {
131  if (m_debug) std::cout << " TrajectoryMeasurement hit Det: Tracker" << std::endl;
132  if (m_debug) std::cout << " TrajectoryMeasurement hit dimension: " << trajMeasurementHitDim << std::endl;
134  double xresid = tsos.localPosition().x() - trajMeasurementHit->localPosition().x();
135  double xresiderr2 = tsos.localError().positionError().xx() + trajMeasurementHit->localPositionError().xx();
136  m_tracker_chi2 += xresid * xresid / xresiderr2;
137 
138 
139 
140  if ( trajMeasurementHitId.subdetId() == StripSubdetector::TID
141  || trajMeasurementHitId.subdetId() == StripSubdetector::TEC) {
142  m_contains_TIDTEC = true;
143  }
144  // YP I add false here. No trajectory measurments in Muon system if we corrected TrackTransformer accordingly
145  } else if ( false && trajMeasurementHitId.det() == DetId::Muon ) {
146 
147  //AR: I removed the false criteria for cosmic tests
148  // } else if (trajMeasurementHitId.det() == DetId::Muon ) {
149 
150  if (m_debug) std::cout << " TrajectoryMeasurement hit Det: Muon" << std::endl;
151 
152  if ( trajMeasurementHitId.subdetId() == MuonSubdetId::DT ) {
153  const DTChamberId chamberId(trajMeasurementHitId.rawId());
154  if (m_debug) std::cout << " TrajectoryMeasurement hit subDet: DT wheel " << chamberId.wheel() << " station " << chamberId.station() << " sector " << chamberId.sector() << std::endl;
155 
156  // double gChX = globalGeometry->idToDet(chamberId)->position().x();
157  // double gChY = globalGeometry->idToDet(chamberId)->position().y();
158  // double gChZ = globalGeometry->idToDet(chamberId)->position().z();
159  // std::cout << " The chamber position in global frame x: " << gChX << " y: " << gChY << " z: " << gChZ << std::endl;
160 
161  const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(chamberId);
162  double chamber_width = geomDet->surface().bounds().width();
163  double chamber_length = geomDet->surface().bounds().length();
164 
165 
166  double hitX = trajMeasurementHit->localPosition().x();
167  double hitY = trajMeasurementHit->localPosition().y();
168  // double hitZ = trajMeasurementHit->localPosition().z();
169 
170  double tsosX = tsos.localPosition().x();
171  double tsosY = tsos.localPosition().y();
172  // double tsosZ = tsos.localPosition().z();
173 
174  int residualDT13IsAdded = false;
175  int residualDT2IsAdded = false;
176 
177 
178  // have we seen this chamber before?
179  if (m_dt13.find(chamberId) == m_dt13.end() && m_dt2.find(chamberId) == m_dt2.end()) {
180  if (m_debug) std::cout << "AR: pushing back chamber: " << chamberId << std::endl;
181  m_chamberIds.push_back(chamberId);
182  //addTrkCovMatrix(chamberId, tsos); // only for the 1st hit
183  }
184  if (m_debug) std::cout << "AR: size of chamberId: " << m_chamberIds.size() << std::endl;
185 
186  if (m_debug) std::cout << " TrajectoryMeasurement hit dimension: " << trajMeasurementHitDim << std::endl;
187  if ( trajMeasurementHitDim > 1 ) {
188  std::vector<const TrackingRecHit*> vDTSeg2D = trajMeasurementHit->recHits();
189  if (m_debug) std::cout << " vDTSeg2D size: " << vDTSeg2D.size() << std::endl;
190  for ( std::vector<const TrackingRecHit*>::const_iterator itDTSeg2D = vDTSeg2D.begin();
191  itDTSeg2D != vDTSeg2D.end();
192  ++itDTSeg2D ) {
193  std::vector<const TrackingRecHit*> vDTHits1D = (*itDTSeg2D)->recHits();
194  if (m_debug) std::cout << " vDTHits1D size: " << vDTHits1D.size() << std::endl;
195  for ( std::vector<const TrackingRecHit*>::const_iterator itDTHits1D = vDTHits1D.begin();
196  itDTHits1D != vDTHits1D.end();
197  ++itDTHits1D ) {
198  const TrackingRecHit* hit = *itDTHits1D;
199  if (m_debug) std::cout << " hit dimension: " << hit->dimension() << std::endl;
200 
201  DetId hitId = hit->geographicalId();
202  const DTSuperLayerId superLayerId(hitId.rawId());
203  const DTLayerId layerId(hitId.rawId());
204  if (m_debug) std::cout << " hit superLayerId: " << superLayerId.superLayer() << std::endl;
205  if (m_debug) std::cout << " hit layerId: " << layerId.layer() << std::endl;
206 
207  if ( superLayerId.superlayer() == 2 && vDTHits1D.size() >= 3 ) {
208  if ( m_dt2.find(chamberId) == m_dt2.end() ) {
209  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
210  m_dt2[chamberId] = new MuonDT2ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
211  if (m_debug) std::cout << " This is first appearance of the DT with hits in superlayer 2" << std::endl;
212 
213  }
214  m_dt2[chamberId]->addResidual(prop, &tsos, hit,chamber_width,chamber_length);
215  residualDT2IsAdded = true;
216 
217  } else if ( (superLayerId.superlayer() == 1 || superLayerId.superlayer() == 3) && vDTHits1D.size() >= 6 ) {
218  if ( m_dt13.find(chamberId) == m_dt13.end() ) {
219  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
220  m_dt13[chamberId] = new MuonDT13ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
221  if (m_debug) std::cout << " This is first appearance of the DT with hits in superlayers 1 and 3" << std::endl;
222  }
223  m_dt13[chamberId]->addResidual(prop, &tsos, hit,chamber_width,chamber_length);
224  residualDT13IsAdded = true;
225 
226 
227  }
228  }
229  }
230  }
231 
232  if ( residualDT13IsAdded ==true && residualDT2IsAdded == true && chamberId.wheel() == 0 && chamberId.station() == 2 && chamberId.sector() == 7 ) {
233  if (m_debug) {
234  std::cout << "MYMARK " << tsosX << " " << hitX << " " << tsosX - hitX << " " << m_dt13[chamberId]->trackx() << " " << m_dt13[chamberId]->residual()
235  << " " << tsosY << " " << hitY << " " << tsosY - hitY << " " << m_dt2[chamberId]->tracky() << " " << m_dt2[chamberId]->residual()
236  << " " << tsosF.localPosition().x() << " " << tsosF.localPosition().y() << " " << tsosF.localPosition().z()
237  << " " << tsosB.localPosition().x() << " " << tsosB.localPosition().y() << " " << tsosB.localPosition().z()
238  << " " << tsosU.localPosition().x() << " " << tsosU.localPosition().y() << " " << tsosU.localPosition().z() << std::endl;
239  }
240  }
241 
242  // http://cmslxr.fnal.gov/lxr/source/DataFormats/TrackReco/src/HitPattern.cc#101
243  // YP I add false here. No trajectory measurments in Muon system if we corrected TrackTransformer accordingly
244  } else if ( false && trajMeasurementHitId.subdetId() == MuonSubdetId::CSC ) {
245  const CSCDetId cscDetId(trajMeasurementHitId.rawId());
246  const CSCDetId chamberId2(cscDetId.endcap(), cscDetId.station(), cscDetId.ring(), cscDetId.chamber());
247  if (m_debug) std::cout << " TrajectoryMeasurement hit subDet: CSC endcap " << cscDetId.endcap() << " station " << cscDetId.station() << " ring " << cscDetId.ring() << " chamber " << cscDetId.chamber() << std::endl;
248  if (m_debug) std::cout << " TrajectoryMeasurement hit dimension: " << trajMeasurementHitDim << std::endl;
249 
250  if ( trajMeasurementHitDim == 4 ) {
251  std::vector<const TrackingRecHit*> vCSCHits2D = trajMeasurementHit->recHits();
252  if (m_debug) std::cout << " vCSCHits2D size: " << vCSCHits2D.size() << std::endl;
253  if ( vCSCHits2D.size() >= 5 ) {
254  for ( std::vector<const TrackingRecHit*>::const_iterator itCSCHits2D = vCSCHits2D.begin();
255  itCSCHits2D != vCSCHits2D.end();
256  ++itCSCHits2D ) {
257  const TrackingRecHit* cscHit2D = *itCSCHits2D;
258  if (m_debug) std::cout << " cscHit2D dimension: " << cscHit2D->dimension() << std::endl;
259  const TrackingRecHit* hit = cscHit2D;
260  if (m_debug) std::cout << " hit dimension: " << hit->dimension() << std::endl;
261 
262  DetId hitId = hit->geographicalId();
263  const CSCDetId cscDetId(hitId.rawId());
264  if (m_debug) std::cout << " hit layer: " << cscDetId.layer() << std::endl;
265 
266  // not sure why we sometimes get layer == 0
267  if (cscDetId.layer() == 0) continue;
268 
269  // have we seen this chamber before?
270  if ( m_csc.find(chamberId2) == m_csc.end() ) {
271  m_chamberIds.push_back(chamberId2);
272  //addTrkCovMatrix(chamberId, tsos); // only for the 1st hit
273  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId2);
274  m_csc[chamberId2] = new MuonCSCChamberResidual(globalGeometry, navigator, chamberId2, chamberAlignable);
275  if (m_debug) std::cout << " This is first appearance of the CSC with hits QQQ" << std::endl;
276  }
277 
278  m_csc[chamberId2]->addResidual(prop, &tsos, hit,250.0,250.0);
279 
280  }
281  }
282  }
283  } else {
284  if (m_debug) std::cout << " TrajectoryMeasurement hit subDet: UNKNOWN" << std::endl;
285  if (m_debug) std::cout << "AR: trajMeasurementHitId.det(): " << trajMeasurementHitId.subdetId() << std::endl;
286  }
287  } else {
288  if (m_debug) std::cout << " TrajectoryMeasurement hit det: UNKNOWN" << std::endl;
289  if (m_debug) std::cout << "AR: trajMeasurementHitId.det(): " << trajMeasurementHitId.det() << std::endl;
290  if (m_debug) std::cout << "DetId::Tracker: " << DetId::Tracker << std::endl;
291  }
292  }
293  }
294 
295 
296 
297 
298  int iT2 = 0, iM2 = 0;
299  for (trackingRecHit_iterator hit2 = m_recoTrack->recHitsBegin(); hit2 != m_recoTrack->recHitsEnd(); ++hit2) {
300  if((*hit2)->isValid()) {
301  DetId hitId2 = (*hit2)->geographicalId();
302  if ( hitId2.det() == DetId::Tracker ) {
303  iT2++;
304  if (m_debug) std::cout << "Tracker Hit " << iT2 << " is found. We don't calcualte Tsos for it" << std::endl;
305  } else if ( hitId2.det() == DetId::Muon ){
306  // if ( (*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit ) {
307  // LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged";
308  // continue;
309  // }
310  iM2++;
311  if (m_debug) std::cout << "Muon Hit " << iM2 << " is found. Dimension: " << (*hit2)->dimension() << std::endl;
312  if ( hitId2.subdetId() == MuonSubdetId::DT ) {
313  const DTChamberId chamberId(hitId2.rawId());
314  if (m_debug) std::cout << "Muon Hit in DT wheel " << chamberId.wheel() << " station " << chamberId.station() << " sector " << chamberId.sector() << std::endl;
315 
316  const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(chamberId);
317  double chamber_width = geomDet->surface().bounds().width();
318  double chamber_length = geomDet->surface().bounds().length();
319 
320 
321 
322  if ( (*hit2)->dimension() > 1 ) {
323  // std::vector<const TrackingRecHit*> vDTSeg2D = (*hit2)->recHits();
324  std::vector<TrackingRecHit*> vDTSeg2D = (*hit2)->recHits();
325 
326  if (m_debug) std::cout << " vDTSeg2D size: " << vDTSeg2D.size() << std::endl;
327 
328  // for ( std::vector<const TrackingRecHit*>::const_iterator itDTSeg2D = vDTSeg2D.begin();
329  // itDTSeg2D != vDTSeg2D.end();
330  // ++itDTSeg2D ) {
331 
332  for ( std::vector<TrackingRecHit*>::const_iterator itDTSeg2D = vDTSeg2D.begin();
333  itDTSeg2D != vDTSeg2D.end();
334  ++itDTSeg2D ) {
335 
336 
337  // std::vector<const TrackingRecHit*> vDTHits1D = (*itDTSeg2D)->recHits();
338  std::vector<TrackingRecHit*> vDTHits1D = (*itDTSeg2D)->recHits();
339  if (m_debug) std::cout << " vDTHits1D size: " << vDTHits1D.size() << std::endl;
340  // for ( std::vector<const TrackingRecHit*>::const_iterator itDTHits1D = vDTHits1D.begin();
341  // itDTHits1D != vDTHits1D.end();
342  // ++itDTHits1D ) {
343  for ( std::vector<TrackingRecHit*>::const_iterator itDTHits1D = vDTHits1D.begin();
344  itDTHits1D != vDTHits1D.end();
345  ++itDTHits1D ) {
346  //const TrackingRecHit* hit = *itDTHits1D;
347  TrackingRecHit* hit = *itDTHits1D;
348  if (m_debug) std::cout << " hit dimension: " << hit->dimension() << std::endl;
349 
350  DetId hitId = hit->geographicalId();
351  const DTSuperLayerId superLayerId(hitId.rawId());
352  const DTLayerId layerId(hitId.rawId());
353  if (m_debug) std::cout << " hit superLayerId: " << superLayerId.superLayer() << std::endl;
354  if (m_debug) std::cout << " hit layerId: " << layerId.layer() << std::endl;
355 
356  if ( superLayerId.superlayer() == 2 && vDTHits1D.size() >= 3 ) {
357  if ( m_dt2.find(chamberId) == m_dt2.end() ) {
358  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
359  m_dt2[chamberId] = new MuonDT2ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
360  if (m_debug) std::cout << " This is first appearance of the DT with hits in superlayer 2" << std::endl;
361 
362  // have we seen this chamber before? check if it was in dt13
363  if ( m_dt13.find(chamberId) == m_dt13.end() ) {
364  m_chamberIds.push_back(chamberId);
365  }
366 
367  }
368 
369  TrajectoryStateOnSurface extrapolation;
370  extrapolation = prop->propagate( lastTrackerTsos, globalGeometry->idToDet(hitId)->surface() );
371 
372  if ( extrapolation.isValid() ) {
373  if (m_debug) {
374  std::cout << " extrapolation localPosition()"
375  << " x: " << extrapolation.localPosition().x()
376  << " y: " << extrapolation.localPosition().y()
377  << " z: " << extrapolation.localPosition().z() << std::endl;
378  }
379  m_dt2[chamberId]->addResidual(prop, &extrapolation, hit,chamber_width,chamber_length);
380  }
381  // residualDT2IsAdded = true;
382 
383  } else if ( (superLayerId.superlayer() == 1 || superLayerId.superlayer() == 3) && vDTHits1D.size() >= 6 ) {
384  if ( m_dt13.find(chamberId) == m_dt13.end() ) {
385  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
386  m_dt13[chamberId] = new MuonDT13ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
387  if (m_debug) std::cout << " This is first appearance of the DT with hits in superlayers 1 and 3" << std::endl;
388 
389  // have we seen this chamber before? check if it was in dt2
390  if ( m_dt2.find(chamberId) == m_dt2.end() ) {
391  m_chamberIds.push_back(chamberId);
392  }
393  }
394 
395  TrajectoryStateOnSurface extrapolation;
396  extrapolation = prop->propagate( lastTrackerTsos, globalGeometry->idToDet(hitId)->surface() );
397 
398  if ( extrapolation.isValid() ) {
399  if (m_debug) {
400  std::cout << " extrapolation localPosition()"
401  << " x: " << extrapolation.localPosition().x()
402  << " y: " << extrapolation.localPosition().y()
403  << " z: " << extrapolation.localPosition().z() << std::endl;
404  }
405  m_dt13[chamberId]->addResidual(prop, &extrapolation, hit,chamber_width,chamber_length);
406  }
407  // residualDT13IsAdded = true;
408 
409 
410  }
411  }
412  }
413  }
414 
415  // std::cout << "Extrapolate last Tracker TSOS to muon hit" << std::endl;
416  // TrajectoryStateOnSurface extrapolation;
417  // extrapolation = prop->propagate( lastTrackerTsos, globalGeometry->idToDet(hitId2)->surface() );
418  //
419  // if ( chamberId2.wheel() == 0 && chamberId2.station() == 2 && chamberId2.sector() == 7 ) {
420  //
421  // double hitX2 = (*hit2)->localPosition().x();
422  // double hitY2 = (*hit2)->localPosition().y();
423  // double hitZ2 = (*hit2)->localPosition().z();
424  //
425  // double tsosX2 = extrapolation.localPosition().x();
426  // double tsosY2 = extrapolation.localPosition().y();
427  // double tsosZ2 = extrapolation.localPosition().z();
428  //
429  // std::cout << "MYMARK " << tsosX2 << " " << hitX2 << " " << tsosX2 - hitX2 << " " << "0" << " " << "0"
430  // << " " << tsosY2 << " " << hitY2 << " " << tsosY2 - hitY2 << " " << "0" << " " << "0"
431  // << " 0 0 0 0 0 0 0 0 0 " << std::endl;
435  // }
436 
437  } else if ( hitId2.subdetId() == MuonSubdetId::CSC ) {
438  const CSCDetId cscDetId2(hitId2.rawId());
439  const CSCDetId chamberId(cscDetId2.endcap(), cscDetId2.station(), cscDetId2.ring(), cscDetId2.chamber());
440  if (m_debug) std::cout << "Muon hit in CSC endcap " << cscDetId2.endcap() << " station " << cscDetId2.station() << " ring " << cscDetId2.ring() << " chamber " << cscDetId2.chamber() << "." << std::endl;
441 
442 
443  if ( (*hit2)->dimension() == 4 ) {
444  // std::vector<const TrackingRecHit*> vCSCHits2D = (*hit2)->recHits();
445  std::vector<TrackingRecHit*> vCSCHits2D = (*hit2)->recHits();
446  if (m_debug) std::cout << " vCSCHits2D size: " << vCSCHits2D.size() << std::endl;
447  if ( vCSCHits2D.size() >= 5 ) {
448  // for ( std::vector<const TrackingRecHit*>::const_iterator itCSCHits2D = vCSCHits2D.begin();
449  // itCSCHits2D != vCSCHits2D.end();
450  // ++itCSCHits2D ) {
451 
452  for ( std::vector<TrackingRecHit*>::const_iterator itCSCHits2D = vCSCHits2D.begin();
453  itCSCHits2D != vCSCHits2D.end();
454  ++itCSCHits2D ) {
455  // const TrackingRecHit* cscHit2D = *itCSCHits2D;
456  TrackingRecHit* cscHit2D = *itCSCHits2D;
457  if (m_debug) std::cout << " cscHit2D dimension: " << cscHit2D->dimension() << std::endl;
458  // const TrackingRecHit* hit = cscHit2D;
459  TrackingRecHit* hit = cscHit2D;
460  if (m_debug) std::cout << " hit dimension: " << hit->dimension() << std::endl;
461 
462  DetId hitId = hit->geographicalId();
463  const CSCDetId cscDetId(hitId.rawId());
464 
465  if (m_debug) {
466  std::cout << " hit layer: " << cscDetId.layer() << std::endl;
467 
468  std::cout << " hit localPosition"
469  << " x: " << hit->localPosition().x()
470  << " y: " << hit->localPosition().y()
471  << " z: " << hit->localPosition().z()
472  << std::endl;
473  std::cout << " hit globalPosition"
474  << " x: " << globalGeometry->idToDet(hitId)->toGlobal(hit->localPosition()).x()
475  << " y: " << globalGeometry->idToDet(hitId)->toGlobal(hit->localPosition()).y()
476  << " z: " << globalGeometry->idToDet(hitId)->toGlobal(hit->localPosition()).z()
477  << std::endl;
478  }
479 
480  // not sure why we sometimes get layer == 0
481  if (cscDetId.layer() == 0) continue;
482 
483  // have we seen this chamber before?
484  if (m_debug) std::cout << "Have we seen this chamber before?";
485  if ( m_csc.find(chamberId) == m_csc.end() ) {
486  if (m_debug) std::cout << " NO. m_csc.count() = " << m_csc.count(chamberId) << std::endl;
487  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
488  m_csc[chamberId] = new MuonCSCChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
489  if (m_debug) std::cout << " This is first appearance of the CSC with hits m_csc.count() = " << m_csc.count(chamberId) << std::endl;
490  m_chamberIds.push_back(chamberId);
491  //addTrkCovMatrix(chamberId, tsos); // only for the 1st hit
492  } else {
493  if (m_debug) std::cout << " YES. m_csc.count() = " << m_csc.count(chamberId) << std::endl;
494  }
495 
496  if (m_debug) {
497  std::cout << " lastTrackerTsos localPosition"
498  << " x: " << lastTrackerTsos.localPosition().x()
499  << " y: " << lastTrackerTsos.localPosition().y()
500  << " z: " << lastTrackerTsos.localPosition().z()
501  << std::endl;
502  std::cout << " lastTrackerTsos globalPosition"
503  << " x: " << lastTrackerTsos.globalPosition().x()
504  << " y: " << lastTrackerTsos.globalPosition().y()
505  << " z: " << lastTrackerTsos.globalPosition().z()
506  << std::endl;
507  std::cout << " Do extrapolation from lastTrackerTsos to hit surface" << std::endl;
508  }
509  TrajectoryStateOnSurface extrapolation;
510  extrapolation = prop->propagate( lastTrackerTsos, globalGeometry->idToDet(hitId)->surface() );
511  if (m_debug) std::cout << " extrapolation.isValid() = " << extrapolation.isValid() << std::endl;
512 
513  if ( extrapolation.isValid() ) {
514  if (m_debug) {
515  std::cout << " extrapolation localPosition()"
516  << " x: " << extrapolation.localPosition().x()
517  << " y: " << extrapolation.localPosition().y()
518  << " z: " << extrapolation.localPosition().z() << std::endl;
519  }
520  m_csc[chamberId]->addResidual(prop, &extrapolation, hit,250.0,250.0);
521  }
522  }
523  }
524  }
525 
526 
527  } else if ( hitId2.subdetId() == MuonSubdetId::RPC ) {
528  if (m_debug) std::cout << "Muon Hit in RPC" << std::endl;
529  } else {
530  if (m_debug) std::cout << "Warning! Muon Hit not in DT or CSC or RPC" << std::endl;
531  }
532  // recHitsForRefit.push_back(theMuonRecHitBuilder->build(&**hit));
533  if ( hitId2.subdetId() == MuonSubdetId::DT || hitId2.subdetId() == MuonSubdetId::CSC ) {
534 
535  }
536  }
537  }
538  }
539 
540 
541 
542  if (m_debug) std::cout << "END MuonResidualsFromTrack" << std::endl << std::endl;
543  }
544 
545 
547  const reco::Muon *recoMuon,
549  double maxResidual )
550  : m_recoMuon(recoMuon)
551  {
552  bool m_debug = false;
553 
554  clear();
557 
558  m_tracker_chi2 = m_recoMuon->innerTrack()->chi2();
559  m_tracker_numHits = m_recoMuon->innerTrack()->ndof() + 5;
561 
562  /*
563  for (trackingRecHit_iterator hit = m_recoMuon->innerTrack()->recHitsBegin(); hit != m_recoMuon->innerTrack()->recHitsEnd(); ++hit)
564  {
565  DetId id = (*hit)->geographicalId();
566  if (id.det() == DetId::Tracker)
567  {
568  m_tracker_numHits++;
569  if (id.subdetId() == StripSubdetector::TID || id.subdetId() == StripSubdetector::TEC) m_contains_TIDTEC = true;
570  }
571  }
572  */
573 
574  for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = m_recoMuon->matches().begin();
575  chamberMatch != m_recoMuon->matches().end(); chamberMatch++)
576  {
577  if (chamberMatch->id.det() != DetId::Muon ) continue;
578 
579  for (std::vector<reco::MuonSegmentMatch>::const_iterator segMatch = chamberMatch->segmentMatches.begin();
580  segMatch != chamberMatch->segmentMatches.end(); ++segMatch)
581  {
582  // select the only segment that belongs to track and is the best in station by dR
583  if (! (segMatch->isMask(reco::MuonSegmentMatch::BestInStationByDR) &&
584  segMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) ) continue;
585 
586  if (chamberMatch->id.subdetId() == MuonSubdetId::DT)
587  {
588  const DTChamberId chamberId(chamberMatch->id.rawId());
589 
590  DTRecSegment4DRef segmentDT = segMatch->dtSegmentRef;
591  const DTRecSegment4D* segment = segmentDT.get();
592  if (segment == nullptr) continue;
593 
594  if ( segment->hasPhi() && fabs(chamberMatch->x - segMatch->x) > maxResidual ) continue;
595  if ( segment->hasZed() && fabs(chamberMatch->y - segMatch->y) > maxResidual ) continue;
596 
597  // have we seen this chamber before?
598  if (m_dt13.find(chamberId) == m_dt13.end() && m_dt2.find(chamberId) == m_dt2.end()) {
599  m_chamberIds.push_back(chamberId);
600  }
601 
602  if (segment->hasZed())
603  {
604  if (m_dt2.find(chamberId) == m_dt2.end())
605  {
606  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
607  // YP
608  // m_dt2[chamberId] = new MuonTrackDT2ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
609  }
610  else if (m_debug) std::cout<<"multi segment match to tmuon: dt2 -- should not happen!"<<std::endl;
611  m_dt2[chamberId]->setSegmentResidual(&(*chamberMatch), &(*segMatch));
612  }
613  if (segment->hasPhi())
614  {
615  if (m_dt13.find(chamberId) == m_dt13.end())
616  {
617  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
618  // YP
619  // m_dt13[chamberId] = new MuonTrackDT13ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
620  }
621  else if (m_debug) std::cout<<"multi segment match to tmuon: dt13 -- should not happen!"<<std::endl;
622  m_dt13[chamberId]->setSegmentResidual(&(*chamberMatch), &(*segMatch));
623  }
624  }
625 
626  else if (chamberMatch->id.subdetId() == MuonSubdetId::CSC)
627  {
628  const CSCDetId cscDetId(chamberMatch->id.rawId());
629  const CSCDetId chamberId(cscDetId.chamberId());
630 
631  if ( fabs(chamberMatch->x - segMatch->x) > maxResidual ) continue;
632 
633  // have we seen this chamber before?
634  if (m_csc.find(chamberId) == m_csc.end())
635  {
636  m_chamberIds.push_back(chamberId);
637  AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
638  // YP
639  // m_csc[chamberId] = new MuonTrackCSCChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
640  }
641  else if (m_debug) std::cout<<"multi segment match to tmuon: csc -- should not happen!"<<std::endl;
642  m_csc[chamberId]->setSegmentResidual(&(*chamberMatch), &(*segMatch));
643  }
644 
645  }
646  }
647  }
648 
649  // This is destructor
650  // It deletes all chambers residulas
652  for (std::map<DetId,MuonChamberResidual*>::const_iterator residual = m_dt13.begin(); residual != m_dt13.end(); ++residual) {
653  delete residual->second;
654  }
655  for (std::map<DetId,MuonChamberResidual*>::const_iterator residual = m_dt2.begin(); residual != m_dt2.end(); ++residual) {
656  delete residual->second;
657  }
658  for (std::map<DetId,MuonChamberResidual*>::const_iterator residual = m_csc.begin(); residual != m_csc.end(); ++residual) {
659  delete residual->second;
660  }
661  }
662 
663 
665  {
666  m_tracker_numHits = 0;
667  m_tracker_chi2 = 0.;
668  m_contains_TIDTEC = false;
669  m_chamberIds.clear();
670  m_dt13.clear();
671  m_dt2.clear();
672  m_csc.clear();
673  m_trkCovMatrix.clear();
674  }
675 
676 
678  {
679  if (m_tracker_numHits > 5) return m_tracker_chi2 / double(m_tracker_numHits - 5);
680  else return -1.;
681  }
682 
683 
685  {
686  if (m_recoMuon) return m_recoTrack->normalizedChi2();
687  return trackerRedChi2();
688  }
689 
690 
692  {
693  if (type == MuonChamberResidual::kDT13) {
694  if (m_dt13.find(chamberId) == m_dt13.end()) return nullptr;
695  return m_dt13[chamberId];
696  }
697  else if (type == MuonChamberResidual::kDT2) {
698  if (m_dt2.find(chamberId) == m_dt2.end()) return nullptr;
699  return m_dt2[chamberId];
700  }
701  else if (type == MuonChamberResidual::kCSC) {
702  if (m_csc.find(chamberId) == m_csc.end()) return nullptr;
703  return m_csc[chamberId];
704  }
705  else return nullptr;
706  }
707 
709  {
710  const AlgebraicSymMatrix55 cov55 = tsos.localError().matrix();
711  TMatrixDSym cov44(4);
712  // change indices from q/p,dxdz,dydz,x,y to x,y,dxdz,dydz
713  int subs[4] = { 3, 4, 1, 2 };
714  for (int i=0;i<4;i++) for (int j=0;j<4;j++) cov44(i,j) = cov55( subs[i], subs[j] );
715  m_trkCovMatrix[chamberId] = cov44;
716  }
717 
719  {
720  bool m_debug = false;
721 
722  TMatrixDSym result(4);
723  if (m_debug) std::cout<<"MuonResidualsFromTrack:: cov initial:"<<std::endl;
724  result.Print();
725  if (m_trkCovMatrix.find(chamberId) == m_trkCovMatrix.end())
726  {
727  if (m_debug) std::cout<<"MuonResidualsFromTrack:: cov does not exist!"<<std::endl;
728  return result;
729  }
730  result = m_trkCovMatrix[chamberId];
731 
732  if (m_debug) std::cout<<"MuonResidualsFromTrack:: cov before:"<<std::endl;
733  result.Print();
734 
735  // add segment's errors in quadratures to track's covariance matrix
736  double r_err;
737  if (m_csc.find(chamberId) == m_csc.end())
738  {
739  r_err = m_csc[chamberId]->residual_error();
740  result(0,0) += r_err*r_err;
741  r_err = m_csc[chamberId]->resslope_error();
742  result(2,2) += r_err*r_err;
743  }
744  if (m_dt13.find(chamberId) == m_dt13.end())
745  {
746  r_err = m_dt13[chamberId]->residual_error();
747  result(0,0) += r_err*r_err;
748  r_err = m_dt13[chamberId]->resslope_error();
749  result(2,2) += r_err*r_err;
750  }
751  if (m_dt2.find(chamberId) == m_dt2.end())
752  {
753  r_err = m_dt2[chamberId]->residual_error();
754  result(1,1) += r_err*r_err;
755  r_err = m_dt2[chamberId]->resslope_error();
756  result(3,3) += r_err*r_err;
757  }
758  if (m_debug) std::cout<<"MuonResidualsFromTrack:: cov after:"<<std::endl;
759  result.Print();
760 
761  return result;
762  }
763 
765  {
766  bool m_debug = false;
767 
768  TMatrixDSym result(4);
769  TMatrixDSym cov44 = covMatrix(chamberId);
770 
771  // invert it using cholesky decomposition
772  TDecompChol decomp(cov44);
773  bool ok = decomp.Invert(result);
774  if (m_debug) std::cout<<"MuonResidualsFromTrack:: corr after:"<<std::endl;
775  result.Print();
776 
777  if (!ok && m_debug) std::cout<<"MuonResidualsFromTrack:: cov inversion failed!"<<std::endl;
778  return result;
779  }
780 
782  {
783  bool m_debug = false;
784 
785  TMatrixD result(4,4);
786  TMatrixDSym corr44 = corrMatrix(chamberId);
787 
788  // get an upper triangular matrix U such that corr = U^T * U
789  TDecompChol decomp(corr44);
790  bool ok = decomp.Decompose();
791  result = decomp.GetU();
792 
793  if (m_debug) std::cout<<"MuonResidualsFromTrack:: corr cholesky after:"<<std::endl;
794  result.Print();
795 
796  if (!ok && m_debug) std::cout<<"MuonResidualsFromTrack:: corr decomposition failed!"<<std::endl;
797  return result;
798  }
type
Definition: HCALResponse.h:21
float xx() const
Definition: LocalError.h:24
virtual float length() const =0
std::vector< DetId > m_chamberIds
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:253
ConstRecHitPointer const & recHit() const
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:561
virtual TrackRef innerTrack() const
Definition: Muon.h:48
const std::string metname
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
std::map< DetId, MuonChamberResidual * > m_dt2
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:47
TMatrixDSym covMatrix(DetId chamberId)
T y() const
Definition: PV3DBase.h:63
MuonResidualsFromTrack(const edm::EventSetup &iSetup, edm::ESHandle< MagneticField > magneticField, edm::ESHandle< GlobalTrackingGeometry > globalGeometry, edm::ESHandle< DetIdAssociator > muonDetIdAssociator_, edm::ESHandle< Propagator > prop, const Trajectory *traj, const reco::Track *recoTrack, AlignableNavigator *navigator, double maxResidual)
GlobalPoint globalPosition() const
std::map< DetId, MuonChamberResidual * > m_csc
const Bounds & bounds() const
Definition: Surface.h:120
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
LocalError positionError() const
static const unsigned int BestInStationByDR
bool isTrackerMuon() const override
Definition: Muon.h:277
virtual float width() const =0
virtual std::vector< const TrackingRecHit * > recHits() const =0
Access to component RecHits (if any)
DataContainer const & measurements() const
Definition: Trajectory.h:196
T sqrt(T t)
Definition: SSEVec.h:18
virtual int dimension() const =0
T z() const
Definition: PV3DBase.h:64
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:41
TMatrixDSym corrMatrix(DetId chamberId)
const AlgebraicSymMatrix55 & matrix() const
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:245
bool hasPhi() const
Does it have the Phi projection?
const LocalTrajectoryError & localError() const
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:104
virtual LocalPoint localPosition() const =0
#define LogTrace(id)
virtual RecHitPointer build(const TrackingRecHit *p) const =0
build a tracking rechit from an existing rechit
TrajectoryStateOnSurface const & forwardPredictedState() const
Access to forward predicted state (from fitter or builder)
std::vector< ConstRecHitPointer > ConstRecHitContainer
virtual const GeomDet * getGeomDet(const DetId &) const =0
bool hasZed() const
Does it have the Z projection?
Definition: DetId.h:18
const reco::Track * m_recoTrack
std::map< DetId, TMatrixDSym > m_trkCovMatrix
TMatrixD choleskyCorrMatrix(DetId chamberId)
static constexpr int RPC
Definition: MuonSubdetId.h:14
bool isValid() const
std::vector< MuonChamberMatch > & matches()
get muon matching information
Definition: Muon.h:144
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
void addTrkCovMatrix(DetId, TrajectoryStateOnSurface &)
static const unsigned int BelongsToTrackByDR
T get() const
Definition: EventSetup.h:63
const GeomDet * idToDet(DetId) const override
TrajectoryStateOnSurface const & updatedState() const
DetId geographicalId() const
static constexpr int DT
Definition: MuonSubdetId.h:12
virtual LocalError localPositionError() const =0
static constexpr int CSC
Definition: MuonSubdetId.h:13
T x() const
Definition: PV3DBase.h:62
std::map< DetId, MuonChamberResidual * > m_dt13
AlignableDetOrUnitPtr alignableFromDetId(const DetId &detid)
Returns AlignableDetOrUnitPtr corresponding to given DetId.
TrajectoryStateOnSurface const & backwardPredictedState() const
Access to backward predicted state (from smoother)
TrajectoryStateCombiner m_tsoscomb
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:39
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:109
MuonChamberResidual * chamberResidual(DetId chamberId, int type)