19 #include "TDecompChol.h" 32 : m_recoTrack(recoTrack)
37 std::cout <<
"BEGIN MuonResidualsFromTrack" << std::endl;
39 LogTrace(metname) <<
"Tracking Component changed!";
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 ==
nullptr)
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;
694 if (
m_dt13.find(chamberId) ==
m_dt13.end())
return nullptr;
698 if (
m_dt2.find(chamberId) ==
m_dt2.end())
return nullptr;
699 return m_dt2[chamberId];
702 if (
m_csc.find(chamberId) ==
m_csc.end())
return nullptr;
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
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
constexpr uint32_t rawId() const
get the raw id
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
bool isTrackerMuon() const override
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
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
auto recHits() const
Access to reconstructed hits on the track.
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
virtual LocalPoint localPosition() const =0
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
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
MuonChamberResidual * chamberResidual(DetId chamberId, int type)