19 #include "TDecompChol.h" 32 : m_recoTrack(recoTrack)
37 std::cout <<
"BEGIN MuonResidualsFromTrack" << std::endl;
39 LogTrace(metname) <<
"Tracking Component changed!";
50 if((*hit)->isValid()) {
51 DetId hitId = (*hit)->geographicalId();
54 if (m_debug)
std::cout <<
"Tracker Hit " << iT <<
" is found. Add to refit. Dimension: " << (*hit)->dimension() << std::endl;
56 recHitsForRefit.push_back( theTrackerRecHitBuilder->
build(&**
hit) );
63 if (m_debug)
std::cout <<
"Muon Hit " << iM <<
" is found. We do not add muon hits to refit. Dimension: " << (*hit)->dimension() << std::endl;
66 if (m_debug)
std::cout <<
"Muon Hit in DT wheel " << chamberId.wheel() <<
" station " << chamberId.station() <<
" sector " << chamberId.sector() <<
"." << std::endl;
69 if (m_debug)
std::cout <<
"Muon hit in CSC endcap " << cscDetId.endcap() <<
" station " << cscDetId.station() <<
" ring " << cscDetId.ring() <<
" chamber " << cscDetId.chamber() <<
"." << std::endl;
71 if (m_debug)
std::cout <<
"Muon Hit in RPC" << std::endl;
73 if (m_debug)
std::cout <<
"Warning! Muon Hit not in DT or CSC or RPC" << std::endl;
86 double lastTrackerTsosGlobalPositionR = 0.0;
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 ) {
95 if (m_debug)
std::cout <<
" TrajectoryMeasurement #" << nTrajMeasurement << std::endl;
103 if (m_debug)
std::cout <<
" TrajectoryMeasurement TSOS validity: " << tsos.
isValid() << std::endl;
107 std::cout <<
" TrajectoryMeasurement TSOS localPosition" 112 std::cout <<
" TrajectoryMeasurement TSOS globalPosition" 115 <<
" R: " << tsosGlobalPositionR
119 if ( tsosGlobalPositionR > lastTrackerTsosGlobalPositionR ) {
120 lastTrackerTsos = tsos;
121 lastTrackerTsosGlobalPositionR = tsosGlobalPositionR;
126 if (m_debug)
std::cout <<
" TrajectoryMeasurement hit validity: " << trajMeasurementHit->
isValid() << std::endl;
127 if ( trajMeasurementHit->
isValid() ) {
129 int trajMeasurementHitDim = trajMeasurementHit->
dimension();
131 if (m_debug)
std::cout <<
" TrajectoryMeasurement hit Det: Tracker" << std::endl;
132 if (m_debug)
std::cout <<
" TrajectoryMeasurement hit dimension: " << trajMeasurementHitDim << std::endl;
145 }
else if (
false && trajMeasurementHitId.
det() ==
DetId::Muon ) {
150 if (m_debug)
std::cout <<
" TrajectoryMeasurement hit Det: Muon" << std::endl;
154 if (m_debug)
std::cout <<
" TrajectoryMeasurement hit subDet: DT wheel " << chamberId.wheel() <<
" station " << chamberId.station() <<
" sector " << chamberId.sector() << std::endl;
174 int residualDT13IsAdded =
false;
175 int residualDT2IsAdded =
false;
180 if (m_debug)
std::cout <<
"AR: pushing back chamber: " << chamberId << std::endl;
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();
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();
204 if (m_debug)
std::cout <<
" hit superLayerId: " << superLayerId.superLayer() << std::endl;
205 if (m_debug)
std::cout <<
" hit layerId: " << layerId.layer() << std::endl;
207 if ( superLayerId.superlayer() == 2 && vDTHits1D.size() >= 3 ) {
211 if (m_debug)
std::cout <<
" This is first appearance of the DT with hits in superlayer 2" << std::endl;
214 m_dt2[chamberId]->addResidual(prop, &tsos, hit,chamber_width,chamber_length);
215 residualDT2IsAdded =
true;
217 }
else if ( (superLayerId.superlayer() == 1 || superLayerId.superlayer() == 3) && vDTHits1D.size() >= 6 ) {
221 if (m_debug)
std::cout <<
" This is first appearance of the DT with hits in superlayers 1 and 3" << std::endl;
223 m_dt13[chamberId]->addResidual(prop, &tsos, hit,chamber_width,chamber_length);
224 residualDT13IsAdded =
true;
232 if ( residualDT13IsAdded ==
true && residualDT2IsAdded ==
true && chamberId.wheel() == 0 && chamberId.station() == 2 && chamberId.sector() == 7 ) {
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()
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;
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();
258 if (m_debug)
std::cout <<
" cscHit2D dimension: " << cscHit2D->
dimension() << std::endl;
264 if (m_debug)
std::cout <<
" hit layer: " << cscDetId.layer() << std::endl;
267 if (cscDetId.layer() == 0)
continue;
270 if (
m_csc.find(chamberId2) ==
m_csc.end() ) {
275 if (m_debug)
std::cout <<
" This is first appearance of the CSC with hits QQQ" << std::endl;
278 m_csc[chamberId2]->addResidual(prop, &tsos, hit,250.0,250.0);
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;
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;
298 int iT2 = 0, iM2 = 0;
300 if((*hit2)->isValid()) {
301 DetId hitId2 = (*hit2)->geographicalId();
304 if (m_debug)
std::cout <<
"Tracker Hit " << iT2 <<
" is found. We don't calcualte Tsos for it" << std::endl;
311 if (m_debug)
std::cout <<
"Muon Hit " << iM2 <<
" is found. Dimension: " << (*hit2)->dimension() << std::endl;
314 if (m_debug)
std::cout <<
"Muon Hit in DT wheel " << chamberId.wheel() <<
" station " << chamberId.station() <<
" sector " << chamberId.sector() << std::endl;
322 if ( (*hit2)->dimension() > 1 ) {
324 std::vector<TrackingRecHit*> vDTSeg2D = (*hit2)->recHits();
326 if (m_debug)
std::cout <<
" vDTSeg2D size: " << vDTSeg2D.size() << std::endl;
332 for ( std::vector<TrackingRecHit*>::const_iterator itDTSeg2D = vDTSeg2D.begin();
333 itDTSeg2D != vDTSeg2D.end();
338 std::vector<TrackingRecHit*> vDTHits1D = (*itDTSeg2D)->recHits();
339 if (m_debug)
std::cout <<
" vDTHits1D size: " << vDTHits1D.size() << std::endl;
343 for ( std::vector<TrackingRecHit*>::const_iterator itDTHits1D = vDTHits1D.begin();
344 itDTHits1D != vDTHits1D.end();
353 if (m_debug)
std::cout <<
" hit superLayerId: " << superLayerId.superLayer() << std::endl;
354 if (m_debug)
std::cout <<
" hit layerId: " << layerId.layer() << std::endl;
356 if ( superLayerId.superlayer() == 2 && vDTHits1D.size() >= 3 ) {
360 if (m_debug)
std::cout <<
" This is first appearance of the DT with hits in superlayer 2" << std::endl;
372 if ( extrapolation.
isValid() ) {
374 std::cout <<
" extrapolation localPosition()" 379 m_dt2[chamberId]->addResidual(prop, &extrapolation, hit,chamber_width,chamber_length);
383 }
else if ( (superLayerId.superlayer() == 1 || superLayerId.superlayer() == 3) && vDTHits1D.size() >= 6 ) {
387 if (m_debug)
std::cout <<
" This is first appearance of the DT with hits in superlayers 1 and 3" << std::endl;
398 if ( extrapolation.
isValid() ) {
400 std::cout <<
" extrapolation localPosition()" 405 m_dt13[chamberId]->addResidual(prop, &extrapolation, hit,chamber_width,chamber_length);
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;
443 if ( (*hit2)->dimension() == 4 ) {
445 std::vector<TrackingRecHit*> vCSCHits2D = (*hit2)->recHits();
446 if (m_debug)
std::cout <<
" vCSCHits2D size: " << vCSCHits2D.size() << std::endl;
447 if ( vCSCHits2D.size() >= 5 ) {
452 for ( std::vector<TrackingRecHit*>::const_iterator itCSCHits2D = vCSCHits2D.begin();
453 itCSCHits2D != vCSCHits2D.end();
457 if (m_debug)
std::cout <<
" cscHit2D dimension: " << cscHit2D->
dimension() << std::endl;
466 std::cout <<
" hit layer: " << cscDetId.layer() << std::endl;
481 if (cscDetId.layer() == 0)
continue;
484 if (m_debug)
std::cout <<
"Have we seen this chamber before?";
486 if (m_debug)
std::cout <<
" NO. m_csc.count() = " <<
m_csc.count(chamberId) << std::endl;
489 if (m_debug)
std::cout <<
" This is first appearance of the CSC with hits m_csc.count() = " <<
m_csc.count(chamberId) << std::endl;
493 if (m_debug)
std::cout <<
" YES. m_csc.count() = " <<
m_csc.count(chamberId) << std::endl;
497 std::cout <<
" lastTrackerTsos localPosition" 502 std::cout <<
" lastTrackerTsos globalPosition" 507 std::cout <<
" Do extrapolation from lastTrackerTsos to hit surface" << std::endl;
511 if (m_debug)
std::cout <<
" extrapolation.isValid() = " << extrapolation.
isValid() << std::endl;
513 if ( extrapolation.
isValid() ) {
515 std::cout <<
" extrapolation localPosition()" 520 m_csc[chamberId]->addResidual(prop, &extrapolation, hit,250.0,250.0);
528 if (m_debug)
std::cout <<
"Muon Hit in RPC" << std::endl;
530 if (m_debug)
std::cout <<
"Warning! Muon Hit not in DT or CSC or RPC" << std::endl;
542 if (m_debug)
std::cout <<
"END MuonResidualsFromTrack" << std::endl << std::endl;
552 bool m_debug =
false;
574 for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch =
m_recoMuon->
matches().begin();
577 if (chamberMatch->id.det() !=
DetId::Muon )
continue;
579 for (std::vector<reco::MuonSegmentMatch>::const_iterator segMatch = chamberMatch->segmentMatches.begin();
580 segMatch != chamberMatch->segmentMatches.end(); ++segMatch)
588 const DTChamberId chamberId(chamberMatch->id.rawId());
592 if (segment == 0)
continue;
594 if ( segment->
hasPhi() && fabs(chamberMatch->x - segMatch->x) > maxResidual )
continue;
595 if ( segment->
hasZed() && fabs(chamberMatch->y - segMatch->y) > maxResidual )
continue;
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));
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));
628 const CSCDetId cscDetId(chamberMatch->id.rawId());
629 const CSCDetId chamberId(cscDetId.chamberId());
631 if ( fabs(chamberMatch->x - segMatch->x) > maxResidual )
continue;
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));
652 for (std::map<DetId,MuonChamberResidual*>::const_iterator residual =
m_dt13.begin(); residual !=
m_dt13.end(); ++residual) {
653 delete residual->second;
655 for (std::map<DetId,MuonChamberResidual*>::const_iterator residual =
m_dt2.begin(); residual !=
m_dt2.end(); ++residual) {
656 delete residual->second;
658 for (std::map<DetId,MuonChamberResidual*>::const_iterator residual =
m_csc.begin(); residual !=
m_csc.end(); ++residual) {
659 delete residual->second;
699 return m_dt2[chamberId];
703 return m_csc[chamberId];
711 TMatrixDSym cov44(4);
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] );
720 bool m_debug =
false;
723 if (m_debug)
std::cout<<
"MuonResidualsFromTrack:: cov initial:"<<std::endl;
727 if (m_debug)
std::cout<<
"MuonResidualsFromTrack:: cov does not exist!"<<std::endl;
732 if (m_debug)
std::cout<<
"MuonResidualsFromTrack:: cov before:"<<std::endl;
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;
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;
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;
758 if (m_debug)
std::cout<<
"MuonResidualsFromTrack:: cov after:"<<std::endl;
766 bool m_debug =
false;
769 TMatrixDSym cov44 =
covMatrix(chamberId);
772 TDecompChol decomp(cov44);
773 bool ok = decomp.Invert(result);
774 if (m_debug)
std::cout<<
"MuonResidualsFromTrack:: corr after:"<<std::endl;
777 if (!ok && m_debug)
std::cout<<
"MuonResidualsFromTrack:: cov inversion failed!"<<std::endl;
783 bool m_debug =
false;
789 TDecompChol decomp(corr44);
790 bool ok = decomp.Decompose();
791 result = decomp.GetU();
793 if (m_debug)
std::cout<<
"MuonResidualsFromTrack:: corr cholesky after:"<<std::endl;
796 if (!ok && m_debug)
std::cout<<
"MuonResidualsFromTrack:: corr decomposition failed!"<<std::endl;
virtual float length() const =0
std::vector< DetId > m_chamberIds
bool isNonnull() const
Checks for non-null.
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)
virtual TrackRef innerTrack() const
const std::string metname
bool isTrackerMuon() const
LocalPoint localPosition() const
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
std::map< DetId, MuonChamberResidual * > m_dt2
TMatrixDSym covMatrix(DetId chamberId)
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
const reco::Muon * m_recoMuon
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
const Plane & surface() const
The nominal surface of the GeomDet.
LocalError positionError() const
static const unsigned int BestInStationByDR
uint32_t rawId() const
get the raw id
virtual float width() const =0
virtual std::vector< const TrackingRecHit * > recHits() const =0
Access to component RecHits (if any)
DataContainer const & measurements() const
virtual int dimension() const =0
TMatrixDSym corrMatrix(DetId chamberId)
const AlgebraicSymMatrix55 & matrix() const
T const * get() const
Returns C++ pointer to the item.
bool hasPhi() const
Does it have the Phi projection?
const LocalTrajectoryError & localError() const
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
virtual LocalPoint localPosition() const =0
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
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
double normalizedChi2() const
virtual const GeomDet * getGeomDet(const DetId &) const =0
bool hasZed() const
Does it have the Z projection?
~MuonResidualsFromTrack()
const reco::Track * m_recoTrack
std::map< DetId, TMatrixDSym > m_trkCovMatrix
TMatrixD choleskyCorrMatrix(DetId chamberId)
double trackerRedChi2() const
std::vector< MuonChamberMatch > & matches()
get muon matching information
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
void addTrkCovMatrix(DetId, TrajectoryStateOnSurface &)
static const unsigned int BelongsToTrackByDR
const GeomDet * idToDet(DetId) const override
TrajectoryStateOnSurface const & updatedState() const
DetId geographicalId() const
virtual LocalError localPositionError() const =0
Detector det() const
get the detector field from this detid
std::map< DetId, MuonChamberResidual * > m_dt13
AlignableDetOrUnitPtr alignableFromDetId(const DetId &detid)
Returns AlignableDetOrUnitPtr corresponding to given DetId.
TrackingRecHitCollection::base::const_iterator trackingRecHit_iterator
iterator over a vector of reference to TrackingRecHit in the same collection
TrajectoryStateOnSurface const & backwardPredictedState() const
Access to backward predicted state (from smoother)
TrajectoryStateCombiner m_tsoscomb
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
MuonChamberResidual * chamberResidual(DetId chamberId, int type)