Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-23 02:04:37 

0001 /*
0002  * $Id: $ 
0003  */
0004 
0005 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
0006 #include "DataFormats/TrackReco/interface/Track.h"
0007 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0008 
0009 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonDT13ChamberResidual.h"
0010 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonDT2ChamberResidual.h"
0011 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonCSCChamberResidual.h"
0012 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonTrackDT13ChamberResidual.h"
0013 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonTrackDT2ChamberResidual.h"
0014 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonTrackCSCChamberResidual.h"
0015 
0016 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0017 #include "TrackingTools/TrackAssociator/interface/DetIdAssociator.h"
0018 
0019 #include "TDecompChol.h"
0020 #include <cmath>
0021 
0022 edm::ESInputTag MuonResidualsFromTrack::builderESInputTag() { return edm::ESInputTag("", "WithTrackAngle"); }
0023 
0024 MuonResidualsFromTrack::MuonResidualsFromTrack(edm::ESHandle<TransientTrackingRecHitBuilder> trackerRecHitBuilder,
0025                                                edm::ESHandle<MagneticField> magneticField,
0026                                                edm::ESHandle<GlobalTrackingGeometry> globalGeometry,
0027                                                edm::ESHandle<DetIdAssociator> muonDetIdAssociator_,
0028                                                edm::ESHandle<Propagator> prop,
0029                                                const Trajectory* traj,
0030                                                const reco::Track* recoTrack,
0031                                                AlignableNavigator* navigator,
0032                                                double maxResidual)
0033     : m_recoTrack(recoTrack) {
0034   init(trackerRecHitBuilder,
0035        magneticField,
0036        globalGeometry,
0037        muonDetIdAssociator_,
0038        prop,
0039        traj,
0040        recoTrack,
0041        navigator,
0042        maxResidual);
0043 }
0044 
0045 void MuonResidualsFromTrack::init(edm::ESHandle<TransientTrackingRecHitBuilder> trackerRecHitBuilder,
0046                                   edm::ESHandle<MagneticField> magneticField,
0047                                   edm::ESHandle<GlobalTrackingGeometry> globalGeometry,
0048                                   edm::ESHandle<DetIdAssociator> muonDetIdAssociator_,
0049                                   edm::ESHandle<Propagator> prop,
0050                                   const Trajectory* traj,
0051                                   const reco::Track* recoTrack,
0052                                   AlignableNavigator* navigator,
0053                                   double maxResidual) {
0054   bool m_debug = false;
0055 
0056   if (m_debug) {
0057     std::cout << "BEGIN MuonResidualsFromTrack" << std::endl;
0058     const std::string metname = " *** MuonResidualsFromTrack *** ";
0059     LogTrace(metname) << "Tracking Component changed!";
0060   }
0061 
0062   clear();
0063 
0064   reco::TransientTrack track(*m_recoTrack, &*magneticField, globalGeometry);
0065   TransientTrackingRecHit::ConstRecHitContainer recHitsForRefit;
0066   int iT = 0, iM = 0;
0067   for (auto const& hit : m_recoTrack->recHits()) {
0068     if (hit->isValid()) {
0069       DetId hitId = hit->geographicalId();
0070       if (hitId.det() == DetId::Tracker) {
0071         iT++;
0072         if (m_debug)
0073           std::cout << "Tracker Hit " << iT << " is found. Add to refit. Dimension: " << hit->dimension() << std::endl;
0074 
0075         recHitsForRefit.push_back(trackerRecHitBuilder->build(&*hit));
0076       } else if (hitId.det() == DetId::Muon) {
0077         //        if ( hit->geographicalId().subdetId() == 3 && !theRPCInTheFit ) {
0078         //          LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged";
0079         //          continue;
0080         //        }
0081         iM++;
0082         if (m_debug)
0083           std::cout << "Muon Hit " << iM
0084                     << " is found. We do not add muon hits to refit. Dimension: " << hit->dimension() << std::endl;
0085         if (hitId.subdetId() == MuonSubdetId::DT) {
0086           const DTChamberId chamberId(hitId.rawId());
0087           if (m_debug)
0088             std::cout << "Muon Hit in DT wheel " << chamberId.wheel() << " station " << chamberId.station()
0089                       << " sector " << chamberId.sector() << "." << std::endl;
0090         } else if (hitId.subdetId() == MuonSubdetId::CSC) {
0091           const CSCDetId cscDetId(hitId.rawId());
0092           if (m_debug)
0093             std::cout << "Muon hit in CSC endcap " << cscDetId.endcap() << " station " << cscDetId.station() << " ring "
0094                       << cscDetId.ring() << " chamber " << cscDetId.chamber() << "." << std::endl;
0095         } else if (hitId.subdetId() == MuonSubdetId::RPC) {
0096           if (m_debug)
0097             std::cout << "Muon Hit in RPC" << std::endl;
0098         } else {
0099           if (m_debug)
0100             std::cout << "Warning! Muon Hit not in DT or CSC or RPC" << std::endl;
0101         }
0102         //        recHitsForRefit.push_back(theMuonRecHitBuilder->build(&*hit));
0103       }
0104     }
0105   }
0106 
0107   //  TrackTransformer trackTransformer();
0108   //  std::vector<Trajectory> vTrackerTrajectory = trackTransformer.transform(track, recHitsForReFit);
0109   //  std::cout << "Tracker trajectories size " << vTrackerTrajectory.size() << std::endl;
0110 
0111   TrajectoryStateOnSurface lastTrackerTsos;
0112   double lastTrackerTsosGlobalPositionR = 0.0;
0113 
0114   std::vector<TrajectoryMeasurement> vTrajMeasurement = traj->measurements();
0115   if (m_debug)
0116     std::cout << "  Size of vector of TrajectoryMeasurements: " << vTrajMeasurement.size() << std::endl;
0117   int nTrajMeasurement = 0;
0118   for (std::vector<TrajectoryMeasurement>::const_iterator iTrajMeasurement = vTrajMeasurement.begin();
0119        iTrajMeasurement != vTrajMeasurement.end();
0120        ++iTrajMeasurement) {
0121     nTrajMeasurement++;
0122     if (m_debug)
0123       std::cout << "    TrajectoryMeasurement #" << nTrajMeasurement << std::endl;
0124 
0125     const TrajectoryMeasurement& trajMeasurement = *iTrajMeasurement;
0126 
0127     TrajectoryStateOnSurface tsos =
0128         m_tsoscomb(trajMeasurement.forwardPredictedState(), trajMeasurement.backwardPredictedState());
0129     const TrajectoryStateOnSurface& tsosF = trajMeasurement.forwardPredictedState();
0130     const TrajectoryStateOnSurface& tsosB = trajMeasurement.backwardPredictedState();
0131     const TrajectoryStateOnSurface& tsosU = trajMeasurement.updatedState();
0132     if (m_debug)
0133       std::cout << "      TrajectoryMeasurement TSOS validity: " << tsos.isValid() << std::endl;
0134     if (tsos.isValid()) {
0135       double tsosGlobalPositionR = sqrt(tsos.globalPosition().x() * tsos.globalPosition().x() +
0136                                         tsos.globalPosition().y() * tsos.globalPosition().y());
0137       if (m_debug) {
0138         std::cout << "         TrajectoryMeasurement TSOS localPosition"
0139                   << " x: " << tsos.localPosition().x() << " y: " << tsos.localPosition().y()
0140                   << " z: " << tsos.localPosition().z() << std::endl;
0141         std::cout << "         TrajectoryMeasurement TSOS globalPosition"
0142                   << " x: " << tsos.globalPosition().x() << " y: " << tsos.globalPosition().y()
0143                   << " R: " << tsosGlobalPositionR << " z: " << tsos.globalPosition().z() << std::endl;
0144       }
0145       if (tsosGlobalPositionR > lastTrackerTsosGlobalPositionR) {
0146         lastTrackerTsos = tsos;
0147         lastTrackerTsosGlobalPositionR = tsosGlobalPositionR;
0148       }
0149     }
0150 
0151     const TransientTrackingRecHit* trajMeasurementHit = &(*trajMeasurement.recHit());
0152     if (m_debug)
0153       std::cout << "      TrajectoryMeasurement hit validity: " << trajMeasurementHit->isValid() << std::endl;
0154     if (trajMeasurementHit->isValid()) {
0155       DetId trajMeasurementHitId = trajMeasurementHit->geographicalId();
0156       int trajMeasurementHitDim = trajMeasurementHit->dimension();
0157       if (trajMeasurementHitId.det() == DetId::Tracker) {
0158         if (m_debug)
0159           std::cout << "      TrajectoryMeasurement hit Det: Tracker" << std::endl;
0160         if (m_debug)
0161           std::cout << "      TrajectoryMeasurement hit dimension: " << trajMeasurementHitDim << std::endl;
0162         m_tracker_numHits++;
0163         double xresid = tsos.localPosition().x() - trajMeasurementHit->localPosition().x();
0164         double xresiderr2 = tsos.localError().positionError().xx() + trajMeasurementHit->localPositionError().xx();
0165         m_tracker_chi2 += xresid * xresid / xresiderr2;
0166 
0167         if (trajMeasurementHitId.subdetId() == StripSubdetector::TID ||
0168             trajMeasurementHitId.subdetId() == StripSubdetector::TEC) {
0169           m_contains_TIDTEC = true;
0170         }
0171         // YP I add false here. No trajectory measurments in Muon system if we corrected TrackTransformer accordingly
0172       } else if (false && trajMeasurementHitId.det() == DetId::Muon) {
0173         //AR: I removed the false criteria for cosmic tests
0174         // } else if (trajMeasurementHitId.det() == DetId::Muon ) {
0175 
0176         if (m_debug)
0177           std::cout << "      TrajectoryMeasurement hit Det: Muon" << std::endl;
0178 
0179         if (trajMeasurementHitId.subdetId() == MuonSubdetId::DT) {
0180           const DTChamberId chamberId(trajMeasurementHitId.rawId());
0181           if (m_debug)
0182             std::cout << "        TrajectoryMeasurement hit subDet: DT wheel " << chamberId.wheel() << " station "
0183                       << chamberId.station() << " sector " << chamberId.sector() << std::endl;
0184 
0185           //          double gChX = globalGeometry->idToDet(chamberId)->position().x();
0186           //          double gChY = globalGeometry->idToDet(chamberId)->position().y();
0187           //          double gChZ = globalGeometry->idToDet(chamberId)->position().z();
0188           //          std::cout << "           The chamber position in global frame x: " << gChX << " y: " << gChY << " z: " << gChZ << std::endl;
0189 
0190           const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(chamberId);
0191           double chamber_width = geomDet->surface().bounds().width();
0192           double chamber_length = geomDet->surface().bounds().length();
0193 
0194           double hitX = trajMeasurementHit->localPosition().x();
0195           double hitY = trajMeasurementHit->localPosition().y();
0196           // double hitZ = trajMeasurementHit->localPosition().z();
0197 
0198           double tsosX = tsos.localPosition().x();
0199           double tsosY = tsos.localPosition().y();
0200           // double tsosZ = tsos.localPosition().z();
0201 
0202           int residualDT13IsAdded = false;
0203           int residualDT2IsAdded = false;
0204 
0205           // have we seen this chamber before?
0206           if (m_dt13.find(chamberId) == m_dt13.end() && m_dt2.find(chamberId) == m_dt2.end()) {
0207             if (m_debug)
0208               std::cout << "AR: pushing back chamber: " << chamberId << std::endl;
0209             m_chamberIds.push_back(chamberId);
0210             //addTrkCovMatrix(chamberId, tsos); // only for the 1st hit
0211           }
0212           if (m_debug)
0213             std::cout << "AR: size of chamberId: " << m_chamberIds.size() << std::endl;
0214 
0215           if (m_debug)
0216             std::cout << "        TrajectoryMeasurement hit dimension: " << trajMeasurementHitDim << std::endl;
0217           if (trajMeasurementHitDim > 1) {
0218             std::vector<const TrackingRecHit*> vDTSeg2D = trajMeasurementHit->recHits();
0219             if (m_debug)
0220               std::cout << "          vDTSeg2D size: " << vDTSeg2D.size() << std::endl;
0221             for (std::vector<const TrackingRecHit*>::const_iterator itDTSeg2D = vDTSeg2D.begin();
0222                  itDTSeg2D != vDTSeg2D.end();
0223                  ++itDTSeg2D) {
0224               std::vector<const TrackingRecHit*> vDTHits1D = (*itDTSeg2D)->recHits();
0225               if (m_debug)
0226                 std::cout << "            vDTHits1D size: " << vDTHits1D.size() << std::endl;
0227               for (std::vector<const TrackingRecHit*>::const_iterator itDTHits1D = vDTHits1D.begin();
0228                    itDTHits1D != vDTHits1D.end();
0229                    ++itDTHits1D) {
0230                 const TrackingRecHit* hit = *itDTHits1D;
0231                 if (m_debug)
0232                   std::cout << "              hit dimension: " << hit->dimension() << std::endl;
0233 
0234                 DetId hitId = hit->geographicalId();
0235                 const DTSuperLayerId superLayerId(hitId.rawId());
0236                 const DTLayerId layerId(hitId.rawId());
0237                 if (m_debug)
0238                   std::cout << "              hit superLayerId: " << superLayerId.superLayer() << std::endl;
0239                 if (m_debug)
0240                   std::cout << "              hit layerId: " << layerId.layer() << std::endl;
0241 
0242                 if (superLayerId.superlayer() == 2 && vDTHits1D.size() >= 3) {
0243                   if (m_dt2.find(chamberId) == m_dt2.end()) {
0244                     AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0245                     m_dt2[chamberId] =
0246                         new MuonDT2ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0247                     if (m_debug)
0248                       std::cout << "              This is first appearance of the DT with hits in superlayer 2"
0249                                 << std::endl;
0250                   }
0251                   m_dt2[chamberId]->addResidual(prop, &tsos, hit, chamber_width, chamber_length);
0252                   residualDT2IsAdded = true;
0253 
0254                 } else if ((superLayerId.superlayer() == 1 || superLayerId.superlayer() == 3) &&
0255                            vDTHits1D.size() >= 6) {
0256                   if (m_dt13.find(chamberId) == m_dt13.end()) {
0257                     AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0258                     m_dt13[chamberId] =
0259                         new MuonDT13ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0260                     if (m_debug)
0261                       std::cout << "              This is first appearance of the DT with hits in superlayers 1 and 3"
0262                                 << std::endl;
0263                   }
0264                   m_dt13[chamberId]->addResidual(prop, &tsos, hit, chamber_width, chamber_length);
0265                   residualDT13IsAdded = true;
0266                 }
0267               }
0268             }
0269           }
0270 
0271           if (residualDT13IsAdded == true && residualDT2IsAdded == true && chamberId.wheel() == 0 &&
0272               chamberId.station() == 2 && chamberId.sector() == 7) {
0273             if (m_debug) {
0274               std::cout << "MYMARK " << tsosX << " " << hitX << " " << tsosX - hitX << " "
0275                         << m_dt13[chamberId]->trackx() << " " << m_dt13[chamberId]->residual() << " " << tsosY << " "
0276                         << hitY << " " << tsosY - hitY << " " << m_dt2[chamberId]->tracky() << " "
0277                         << m_dt2[chamberId]->residual() << " " << tsosF.localPosition().x() << " "
0278                         << tsosF.localPosition().y() << " " << tsosF.localPosition().z() << " "
0279                         << tsosB.localPosition().x() << " " << tsosB.localPosition().y() << " "
0280                         << tsosB.localPosition().z() << " " << tsosU.localPosition().x() << " "
0281                         << tsosU.localPosition().y() << " " << tsosU.localPosition().z() << std::endl;
0282             }
0283           }
0284 
0285           // http://cmslxr.fnal.gov/lxr/source/DataFormats/TrackReco/src/HitPattern.cc#101
0286           // YP I add false here. No trajectory measurments in Muon system if we corrected TrackTransformer accordingly
0287         } else if (false && trajMeasurementHitId.subdetId() == MuonSubdetId::CSC) {
0288           const CSCDetId cscDetId(trajMeasurementHitId.rawId());
0289           const CSCDetId chamberId2(cscDetId.endcap(), cscDetId.station(), cscDetId.ring(), cscDetId.chamber());
0290           if (m_debug)
0291             std::cout << "        TrajectoryMeasurement hit subDet: CSC endcap " << cscDetId.endcap() << " station "
0292                       << cscDetId.station() << " ring " << cscDetId.ring() << " chamber " << cscDetId.chamber()
0293                       << std::endl;
0294           if (m_debug)
0295             std::cout << "        TrajectoryMeasurement hit dimension: " << trajMeasurementHitDim << std::endl;
0296 
0297           if (trajMeasurementHitDim == 4) {
0298             std::vector<const TrackingRecHit*> vCSCHits2D = trajMeasurementHit->recHits();
0299             if (m_debug)
0300               std::cout << "          vCSCHits2D size: " << vCSCHits2D.size() << std::endl;
0301             if (vCSCHits2D.size() >= 5) {
0302               for (std::vector<const TrackingRecHit*>::const_iterator itCSCHits2D = vCSCHits2D.begin();
0303                    itCSCHits2D != vCSCHits2D.end();
0304                    ++itCSCHits2D) {
0305                 const TrackingRecHit* cscHit2D = *itCSCHits2D;
0306                 if (m_debug)
0307                   std::cout << "            cscHit2D dimension: " << cscHit2D->dimension() << std::endl;
0308                 const TrackingRecHit* hit = cscHit2D;
0309                 if (m_debug)
0310                   std::cout << "              hit dimension: " << hit->dimension() << std::endl;
0311 
0312                 DetId hitId = hit->geographicalId();
0313                 const CSCDetId cscDetId(hitId.rawId());
0314                 if (m_debug)
0315                   std::cout << "              hit layer: " << cscDetId.layer() << std::endl;
0316 
0317                 // not sure why we sometimes get layer == 0
0318                 if (cscDetId.layer() == 0)
0319                   continue;
0320 
0321                 // have we seen this chamber before?
0322                 if (m_csc.find(chamberId2) == m_csc.end()) {
0323                   m_chamberIds.push_back(chamberId2);
0324                   //addTrkCovMatrix(chamberId, tsos); // only for the 1st hit
0325                   AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId2);
0326                   m_csc[chamberId2] =
0327                       new MuonCSCChamberResidual(globalGeometry, navigator, chamberId2, chamberAlignable);
0328                   if (m_debug)
0329                     std::cout << "              This is first appearance of the CSC with hits QQQ" << std::endl;
0330                 }
0331 
0332                 m_csc[chamberId2]->addResidual(prop, &tsos, hit, 250.0, 250.0);
0333               }
0334             }
0335           }
0336         } else {
0337           if (m_debug)
0338             std::cout << "        TrajectoryMeasurement hit subDet: UNKNOWN" << std::endl;
0339           if (m_debug)
0340             std::cout << "AR: trajMeasurementHitId.det(): " << trajMeasurementHitId.subdetId() << std::endl;
0341         }
0342       } else {
0343         if (m_debug)
0344           std::cout << "      TrajectoryMeasurement hit det: UNKNOWN" << std::endl;
0345         if (m_debug)
0346           std::cout << "AR: trajMeasurementHitId.det(): " << trajMeasurementHitId.det() << std::endl;
0347         if (m_debug)
0348           std::cout << "DetId::Tracker: " << DetId::Tracker << std::endl;
0349       }
0350     }
0351   }
0352 
0353   int iT2 = 0, iM2 = 0;
0354   for (auto const& hit2 : m_recoTrack->recHits()) {
0355     if (hit2->isValid()) {
0356       DetId hitId2 = hit2->geographicalId();
0357       if (hitId2.det() == DetId::Tracker) {
0358         iT2++;
0359         if (m_debug)
0360           std::cout << "Tracker Hit " << iT2 << " is found. We don't calcualte Tsos for it" << std::endl;
0361       } else if (hitId2.det() == DetId::Muon) {
0362         //        if ( (*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit ) {
0363         //          LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged";
0364         //          continue;
0365         //        }
0366         iM2++;
0367         if (m_debug)
0368           std::cout << "Muon Hit " << iM2 << " is found. Dimension: " << hit2->dimension() << std::endl;
0369         if (hitId2.subdetId() == MuonSubdetId::DT) {
0370           const DTChamberId chamberId(hitId2.rawId());
0371           if (m_debug)
0372             std::cout << "Muon Hit in DT wheel " << chamberId.wheel() << " station " << chamberId.station()
0373                       << " sector " << chamberId.sector() << std::endl;
0374 
0375           const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(chamberId);
0376           double chamber_width = geomDet->surface().bounds().width();
0377           double chamber_length = geomDet->surface().bounds().length();
0378 
0379           if (hit2->dimension() > 1) {
0380             // std::vector<const TrackingRecHit*> vDTSeg2D = hit2->recHits();
0381             std::vector<TrackingRecHit*> vDTSeg2D = hit2->recHits();
0382 
0383             if (m_debug)
0384               std::cout << "          vDTSeg2D size: " << vDTSeg2D.size() << std::endl;
0385 
0386             // for ( std::vector<const TrackingRecHit*>::const_iterator itDTSeg2D =  vDTSeg2D.begin();
0387             //                                                          itDTSeg2D != vDTSeg2D.end();
0388             //                                                        ++itDTSeg2D ) {
0389 
0390             for (std::vector<TrackingRecHit*>::const_iterator itDTSeg2D = vDTSeg2D.begin(); itDTSeg2D != vDTSeg2D.end();
0391                  ++itDTSeg2D) {
0392               // std::vector<const TrackingRecHit*> vDTHits1D =  (*itDTSeg2D)->recHits();
0393               std::vector<TrackingRecHit*> vDTHits1D = (*itDTSeg2D)->recHits();
0394               if (m_debug)
0395                 std::cout << "            vDTHits1D size: " << vDTHits1D.size() << std::endl;
0396               // for ( std::vector<const TrackingRecHit*>::const_iterator itDTHits1D =  vDTHits1D.begin();
0397               //                                                          itDTHits1D != vDTHits1D.end();
0398               //                                                        ++itDTHits1D ) {
0399               for (std::vector<TrackingRecHit*>::const_iterator itDTHits1D = vDTHits1D.begin();
0400                    itDTHits1D != vDTHits1D.end();
0401                    ++itDTHits1D) {
0402                 //const TrackingRecHit* hit = *itDTHits1D;
0403                 TrackingRecHit* hit = *itDTHits1D;
0404                 if (m_debug)
0405                   std::cout << "              hit dimension: " << hit->dimension() << std::endl;
0406 
0407                 DetId hitId = hit->geographicalId();
0408                 const DTSuperLayerId superLayerId(hitId.rawId());
0409                 const DTLayerId layerId(hitId.rawId());
0410                 if (m_debug)
0411                   std::cout << "              hit superLayerId: " << superLayerId.superLayer() << std::endl;
0412                 if (m_debug)
0413                   std::cout << "              hit layerId: " << layerId.layer() << std::endl;
0414 
0415                 if (superLayerId.superlayer() == 2 && vDTHits1D.size() >= 3) {
0416                   if (m_dt2.find(chamberId) == m_dt2.end()) {
0417                     AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0418                     m_dt2[chamberId] =
0419                         new MuonDT2ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0420                     if (m_debug)
0421                       std::cout << "              This is first appearance of the DT with hits in superlayer 2"
0422                                 << std::endl;
0423 
0424                     // have we seen this chamber before? check if it was in dt13
0425                     if (m_dt13.find(chamberId) == m_dt13.end()) {
0426                       m_chamberIds.push_back(chamberId);
0427                     }
0428                   }
0429 
0430                   TrajectoryStateOnSurface extrapolation;
0431                   extrapolation = prop->propagate(lastTrackerTsos, globalGeometry->idToDet(hitId)->surface());
0432 
0433                   if (extrapolation.isValid()) {
0434                     if (m_debug) {
0435                       std::cout << " extrapolation localPosition()"
0436                                 << " x: " << extrapolation.localPosition().x()
0437                                 << " y: " << extrapolation.localPosition().y()
0438                                 << " z: " << extrapolation.localPosition().z() << std::endl;
0439                     }
0440                     m_dt2[chamberId]->addResidual(prop, &extrapolation, hit, chamber_width, chamber_length);
0441                   }
0442                   //                    residualDT2IsAdded = true;
0443 
0444                 } else if ((superLayerId.superlayer() == 1 || superLayerId.superlayer() == 3) &&
0445                            vDTHits1D.size() >= 6) {
0446                   if (m_dt13.find(chamberId) == m_dt13.end()) {
0447                     AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0448                     m_dt13[chamberId] =
0449                         new MuonDT13ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0450                     if (m_debug)
0451                       std::cout << "              This is first appearance of the DT with hits in superlayers 1 and 3"
0452                                 << std::endl;
0453 
0454                     // have we seen this chamber before? check if it was in dt2
0455                     if (m_dt2.find(chamberId) == m_dt2.end()) {
0456                       m_chamberIds.push_back(chamberId);
0457                     }
0458                   }
0459 
0460                   TrajectoryStateOnSurface extrapolation;
0461                   extrapolation = prop->propagate(lastTrackerTsos, globalGeometry->idToDet(hitId)->surface());
0462 
0463                   if (extrapolation.isValid()) {
0464                     if (m_debug) {
0465                       std::cout << " extrapolation localPosition()"
0466                                 << " x: " << extrapolation.localPosition().x()
0467                                 << " y: " << extrapolation.localPosition().y()
0468                                 << " z: " << extrapolation.localPosition().z() << std::endl;
0469                     }
0470                     m_dt13[chamberId]->addResidual(prop, &extrapolation, hit, chamber_width, chamber_length);
0471                   }
0472                   //                    residualDT13IsAdded = true;
0473                 }
0474               }
0475             }
0476           }
0477 
0478           //          std::cout << "Extrapolate last Tracker TSOS to muon hit" << std::endl;
0479           //            TrajectoryStateOnSurface extrapolation;
0480           //            extrapolation = prop->propagate( lastTrackerTsos, globalGeometry->idToDet(hitId2)->surface() );
0481           //
0482           //            if ( chamberId2.wheel() == 0 && chamberId2.station() == 2 && chamberId2.sector() == 7 ) {
0483           //
0484           //            double hitX2 = hit2->localPosition().x();
0485           //          double hitY2 = hit2->localPosition().y();
0486           //          double hitZ2 = hit2->localPosition().z();
0487           //
0488           //          double tsosX2 = extrapolation.localPosition().x();
0489           //          double tsosY2 = extrapolation.localPosition().y();
0490           //          double tsosZ2 = extrapolation.localPosition().z();
0491           //
0492           //            std::cout << "MYMARK " << tsosX2 << " " << hitX2 << " " << tsosX2 - hitX2 << " " << "0" << " " << "0"
0493           //                            << " " << tsosY2 << " " << hitY2 << " " << tsosY2 - hitY2 << " " << "0"  << " " << "0"
0494           //                            << " 0 0 0 0 0 0 0 0 0 " << std::endl;
0495           ////                            << " " << tsosF.localPosition().x() << " " << tsosF.localPosition().y() << " " << tsosF.localPosition().z()
0496           ////                            << " " << tsosB.localPosition().x() << " " << tsosB.localPosition().y() << " " << tsosB.localPosition().z()
0497           ////                            << " " << tsosU.localPosition().x() << " " << tsosU.localPosition().y() << " " << tsosU.localPosition().z() << std::endl;
0498           //            }
0499 
0500         } else if (hitId2.subdetId() == MuonSubdetId::CSC) {
0501           const CSCDetId cscDetId2(hitId2.rawId());
0502           const CSCDetId chamberId(cscDetId2.endcap(), cscDetId2.station(), cscDetId2.ring(), cscDetId2.chamber());
0503           if (m_debug)
0504             std::cout << "Muon hit in CSC endcap " << cscDetId2.endcap() << " station " << cscDetId2.station()
0505                       << " ring " << cscDetId2.ring() << " chamber " << cscDetId2.chamber() << "." << std::endl;
0506 
0507           if (hit2->dimension() == 4) {
0508             // std::vector<const TrackingRecHit*> vCSCHits2D = hit2->recHits();
0509             std::vector<TrackingRecHit*> vCSCHits2D = hit2->recHits();
0510             if (m_debug)
0511               std::cout << "          vCSCHits2D size: " << vCSCHits2D.size() << std::endl;
0512             if (vCSCHits2D.size() >= 5) {
0513               // for ( std::vector<const TrackingRecHit*>::const_iterator itCSCHits2D =  vCSCHits2D.begin();
0514               //                                                          itCSCHits2D != vCSCHits2D.end();
0515               //                                                        ++itCSCHits2D ) {
0516 
0517               for (std::vector<TrackingRecHit*>::const_iterator itCSCHits2D = vCSCHits2D.begin();
0518                    itCSCHits2D != vCSCHits2D.end();
0519                    ++itCSCHits2D) {
0520                 // const TrackingRecHit* cscHit2D = *itCSCHits2D;
0521                 TrackingRecHit* cscHit2D = *itCSCHits2D;
0522                 if (m_debug)
0523                   std::cout << "            cscHit2D dimension: " << cscHit2D->dimension() << std::endl;
0524                 // const TrackingRecHit* hit = cscHit2D;
0525                 TrackingRecHit* hit = cscHit2D;
0526                 if (m_debug)
0527                   std::cout << "              hit dimension: " << hit->dimension() << std::endl;
0528 
0529                 DetId hitId = hit->geographicalId();
0530                 const CSCDetId cscDetId(hitId.rawId());
0531 
0532                 if (m_debug) {
0533                   std::cout << "              hit layer: " << cscDetId.layer() << std::endl;
0534 
0535                   std::cout << " hit localPosition"
0536                             << " x: " << hit->localPosition().x() << " y: " << hit->localPosition().y()
0537                             << " z: " << hit->localPosition().z() << std::endl;
0538                   std::cout << " hit globalPosition"
0539                             << " x: " << globalGeometry->idToDet(hitId)->toGlobal(hit->localPosition()).x()
0540                             << " y: " << globalGeometry->idToDet(hitId)->toGlobal(hit->localPosition()).y()
0541                             << " z: " << globalGeometry->idToDet(hitId)->toGlobal(hit->localPosition()).z()
0542                             << std::endl;
0543                 }
0544 
0545                 // not sure why we sometimes get layer == 0
0546                 if (cscDetId.layer() == 0)
0547                   continue;
0548 
0549                 // have we seen this chamber before?
0550                 if (m_debug)
0551                   std::cout << "Have we seen this chamber before?";
0552                 if (m_csc.find(chamberId) == m_csc.end()) {
0553                   if (m_debug)
0554                     std::cout << " NO. m_csc.count() = " << m_csc.count(chamberId) << std::endl;
0555                   AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0556                   m_csc[chamberId] = new MuonCSCChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0557                   if (m_debug)
0558                     std::cout << "              This is first appearance of the CSC with hits m_csc.count() = "
0559                               << m_csc.count(chamberId) << std::endl;
0560                   m_chamberIds.push_back(chamberId);
0561                   //addTrkCovMatrix(chamberId, tsos); // only for the 1st hit
0562                 } else {
0563                   if (m_debug)
0564                     std::cout << " YES. m_csc.count() = " << m_csc.count(chamberId) << std::endl;
0565                 }
0566 
0567                 if (m_debug) {
0568                   std::cout << " lastTrackerTsos localPosition"
0569                             << " x: " << lastTrackerTsos.localPosition().x()
0570                             << " y: " << lastTrackerTsos.localPosition().y()
0571                             << " z: " << lastTrackerTsos.localPosition().z() << std::endl;
0572                   std::cout << " lastTrackerTsos globalPosition"
0573                             << " x: " << lastTrackerTsos.globalPosition().x()
0574                             << " y: " << lastTrackerTsos.globalPosition().y()
0575                             << " z: " << lastTrackerTsos.globalPosition().z() << std::endl;
0576                   std::cout << " Do extrapolation from lastTrackerTsos to hit surface" << std::endl;
0577                 }
0578                 TrajectoryStateOnSurface extrapolation;
0579                 extrapolation = prop->propagate(lastTrackerTsos, globalGeometry->idToDet(hitId)->surface());
0580                 if (m_debug)
0581                   std::cout << " extrapolation.isValid() = " << extrapolation.isValid() << std::endl;
0582 
0583                 if (extrapolation.isValid()) {
0584                   if (m_debug) {
0585                     std::cout << " extrapolation localPosition()"
0586                               << " x: " << extrapolation.localPosition().x()
0587                               << " y: " << extrapolation.localPosition().y()
0588                               << " z: " << extrapolation.localPosition().z() << std::endl;
0589                   }
0590                   m_csc[chamberId]->addResidual(prop, &extrapolation, hit, 250.0, 250.0);
0591                 }
0592               }
0593             }
0594           }
0595 
0596         } else if (hitId2.subdetId() == MuonSubdetId::RPC) {
0597           if (m_debug)
0598             std::cout << "Muon Hit in RPC" << std::endl;
0599         } else {
0600           if (m_debug)
0601             std::cout << "Warning! Muon Hit not in DT or CSC or RPC" << std::endl;
0602         }
0603         //        recHitsForRefit.push_back(theMuonRecHitBuilder->build(&**hit));
0604         if (hitId2.subdetId() == MuonSubdetId::DT || hitId2.subdetId() == MuonSubdetId::CSC) {
0605         }
0606       }
0607     }
0608   }
0609 
0610   if (m_debug)
0611     std::cout << "END MuonResidualsFromTrack" << std::endl << std::endl;
0612 }
0613 
0614 MuonResidualsFromTrack::MuonResidualsFromTrack(edm::ESHandle<GlobalTrackingGeometry> globalGeometry,
0615                                                const reco::Muon* recoMuon,
0616                                                AlignableNavigator* navigator,
0617                                                double maxResidual)
0618     : m_recoMuon(recoMuon) {
0619   bool m_debug = false;
0620 
0621   clear();
0622   assert(m_recoMuon->isTrackerMuon() && m_recoMuon->innerTrack().isNonnull());
0623   m_recoTrack = m_recoMuon->innerTrack().get();
0624 
0625   m_tracker_chi2 = m_recoMuon->innerTrack()->chi2();
0626   m_tracker_numHits = m_recoMuon->innerTrack()->ndof() + 5;
0627   m_tracker_numHits = m_tracker_numHits > 0 ? m_tracker_numHits : 0;
0628 
0629   /*
0630                    for(auto const& hit : m_recoMuon->innerTrack()->recHits())
0631                    {
0632                    DetId id = hit->geographicalId();
0633                    if (id.det() == DetId::Tracker)
0634                    {
0635                    m_tracker_numHits++;
0636                    if (id.subdetId() == StripSubdetector::TID  ||  id.subdetId() == StripSubdetector::TEC) m_contains_TIDTEC = true;
0637                    }
0638                    }
0639                    */
0640 
0641   for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = m_recoMuon->matches().begin();
0642        chamberMatch != m_recoMuon->matches().end();
0643        chamberMatch++) {
0644     if (chamberMatch->id.det() != DetId::Muon)
0645       continue;
0646 
0647     for (std::vector<reco::MuonSegmentMatch>::const_iterator segMatch = chamberMatch->segmentMatches.begin();
0648          segMatch != chamberMatch->segmentMatches.end();
0649          ++segMatch) {
0650       // select the only segment that belongs to track and is the best in station by dR
0651       if (!(segMatch->isMask(reco::MuonSegmentMatch::BestInStationByDR) &&
0652             segMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)))
0653         continue;
0654 
0655       if (chamberMatch->id.subdetId() == MuonSubdetId::DT) {
0656         const DTChamberId chamberId(chamberMatch->id.rawId());
0657 
0658         DTRecSegment4DRef segmentDT = segMatch->dtSegmentRef;
0659         const DTRecSegment4D* segment = segmentDT.get();
0660         if (segment == nullptr)
0661           continue;
0662 
0663         if (segment->hasPhi() && fabs(chamberMatch->x - segMatch->x) > maxResidual)
0664           continue;
0665         if (segment->hasZed() && fabs(chamberMatch->y - segMatch->y) > maxResidual)
0666           continue;
0667 
0668         // have we seen this chamber before?
0669         if (m_dt13.find(chamberId) == m_dt13.end() && m_dt2.find(chamberId) == m_dt2.end()) {
0670           m_chamberIds.push_back(chamberId);
0671         }
0672 
0673         if (segment->hasZed()) {
0674           if (m_dt2.find(chamberId) == m_dt2.end()) {
0675             AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0676             // YP
0677             //            m_dt2[chamberId] = new MuonTrackDT2ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0678           } else if (m_debug)
0679             std::cout << "multi segment match to tmuon: dt2  -- should not happen!" << std::endl;
0680           m_dt2[chamberId]->setSegmentResidual(&(*chamberMatch), &(*segMatch));
0681         }
0682         if (segment->hasPhi()) {
0683           if (m_dt13.find(chamberId) == m_dt13.end()) {
0684             AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0685             // YP
0686             //            m_dt13[chamberId] = new MuonTrackDT13ChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0687           } else if (m_debug)
0688             std::cout << "multi segment match to tmuon: dt13  -- should not happen!" << std::endl;
0689           m_dt13[chamberId]->setSegmentResidual(&(*chamberMatch), &(*segMatch));
0690         }
0691       }
0692 
0693       else if (chamberMatch->id.subdetId() == MuonSubdetId::CSC) {
0694         const CSCDetId cscDetId(chamberMatch->id.rawId());
0695         const CSCDetId chamberId(cscDetId.chamberId());
0696 
0697         if (fabs(chamberMatch->x - segMatch->x) > maxResidual)
0698           continue;
0699 
0700         // have we seen this chamber before?
0701         if (m_csc.find(chamberId) == m_csc.end()) {
0702           m_chamberIds.push_back(chamberId);
0703           AlignableDetOrUnitPtr chamberAlignable = navigator->alignableFromDetId(chamberId);
0704           // YP
0705           //          m_csc[chamberId] = new MuonTrackCSCChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable);
0706         } else if (m_debug)
0707           std::cout << "multi segment match to tmuon: csc  -- should not happen!" << std::endl;
0708         m_csc[chamberId]->setSegmentResidual(&(*chamberMatch), &(*segMatch));
0709       }
0710     }
0711   }
0712 }
0713 
0714 // This is destructor
0715 // It deletes all chambers residulas
0716 MuonResidualsFromTrack::~MuonResidualsFromTrack() {
0717   for (std::map<DetId, MuonChamberResidual*>::const_iterator residual = m_dt13.begin(); residual != m_dt13.end();
0718        ++residual) {
0719     delete residual->second;
0720   }
0721   for (std::map<DetId, MuonChamberResidual*>::const_iterator residual = m_dt2.begin(); residual != m_dt2.end();
0722        ++residual) {
0723     delete residual->second;
0724   }
0725   for (std::map<DetId, MuonChamberResidual*>::const_iterator residual = m_csc.begin(); residual != m_csc.end();
0726        ++residual) {
0727     delete residual->second;
0728   }
0729 }
0730 
0731 void MuonResidualsFromTrack::clear() {
0732   m_tracker_numHits = 0;
0733   m_tracker_chi2 = 0.;
0734   m_contains_TIDTEC = false;
0735   m_chamberIds.clear();
0736   m_dt13.clear();
0737   m_dt2.clear();
0738   m_csc.clear();
0739   m_trkCovMatrix.clear();
0740 }
0741 
0742 double MuonResidualsFromTrack::trackerRedChi2() const {
0743   if (m_tracker_numHits > 5)
0744     return m_tracker_chi2 / double(m_tracker_numHits - 5);
0745   else
0746     return -1.;
0747 }
0748 
0749 double MuonResidualsFromTrack::normalizedChi2() const {
0750   if (m_recoMuon)
0751     return m_recoTrack->normalizedChi2();
0752   return trackerRedChi2();
0753 }
0754 
0755 MuonChamberResidual* MuonResidualsFromTrack::chamberResidual(DetId chamberId, int type) {
0756   if (type == MuonChamberResidual::kDT13) {
0757     if (m_dt13.find(chamberId) == m_dt13.end())
0758       return nullptr;
0759     return m_dt13[chamberId];
0760   } else if (type == MuonChamberResidual::kDT2) {
0761     if (m_dt2.find(chamberId) == m_dt2.end())
0762       return nullptr;
0763     return m_dt2[chamberId];
0764   } else if (type == MuonChamberResidual::kCSC) {
0765     if (m_csc.find(chamberId) == m_csc.end())
0766       return nullptr;
0767     return m_csc[chamberId];
0768   } else
0769     return nullptr;
0770 }
0771 
0772 void MuonResidualsFromTrack::addTrkCovMatrix(DetId chamberId, TrajectoryStateOnSurface& tsos) {
0773   const AlgebraicSymMatrix55 cov55 = tsos.localError().matrix();
0774   TMatrixDSym cov44(4);
0775   // change indices from q/p,dxdz,dydz,x,y   to   x,y,dxdz,dydz
0776   int subs[4] = {3, 4, 1, 2};
0777   for (int i = 0; i < 4; i++)
0778     for (int j = 0; j < 4; j++)
0779       cov44(i, j) = cov55(subs[i], subs[j]);
0780   m_trkCovMatrix[chamberId] = cov44;
0781 }
0782 
0783 TMatrixDSym MuonResidualsFromTrack::covMatrix(DetId chamberId) {
0784   bool m_debug = false;
0785 
0786   TMatrixDSym result(4);
0787   if (m_debug)
0788     std::cout << "MuonResidualsFromTrack:: cov initial:" << std::endl;
0789   result.Print();
0790   if (m_trkCovMatrix.find(chamberId) == m_trkCovMatrix.end()) {
0791     if (m_debug)
0792       std::cout << "MuonResidualsFromTrack:: cov does not exist!" << std::endl;
0793     return result;
0794   }
0795   result = m_trkCovMatrix[chamberId];
0796 
0797   if (m_debug)
0798     std::cout << "MuonResidualsFromTrack:: cov before:" << std::endl;
0799   result.Print();
0800 
0801   // add segment's errors in quadratures to track's covariance matrix
0802   double r_err;
0803   if (m_csc.find(chamberId) == m_csc.end()) {
0804     r_err = m_csc[chamberId]->residual_error();
0805     result(0, 0) += r_err * r_err;
0806     r_err = m_csc[chamberId]->resslope_error();
0807     result(2, 2) += r_err * r_err;
0808   }
0809   if (m_dt13.find(chamberId) == m_dt13.end()) {
0810     r_err = m_dt13[chamberId]->residual_error();
0811     result(0, 0) += r_err * r_err;
0812     r_err = m_dt13[chamberId]->resslope_error();
0813     result(2, 2) += r_err * r_err;
0814   }
0815   if (m_dt2.find(chamberId) == m_dt2.end()) {
0816     r_err = m_dt2[chamberId]->residual_error();
0817     result(1, 1) += r_err * r_err;
0818     r_err = m_dt2[chamberId]->resslope_error();
0819     result(3, 3) += r_err * r_err;
0820   }
0821   if (m_debug)
0822     std::cout << "MuonResidualsFromTrack:: cov after:" << std::endl;
0823   result.Print();
0824 
0825   return result;
0826 }
0827 
0828 TMatrixDSym MuonResidualsFromTrack::corrMatrix(DetId chamberId) {
0829   bool m_debug = false;
0830 
0831   TMatrixDSym result(4);
0832   TMatrixDSym cov44 = covMatrix(chamberId);
0833 
0834   // invert it using cholesky decomposition
0835   TDecompChol decomp(cov44);
0836   bool ok = decomp.Invert(result);
0837   if (m_debug)
0838     std::cout << "MuonResidualsFromTrack:: corr after:" << std::endl;
0839   result.Print();
0840 
0841   if (!ok && m_debug)
0842     std::cout << "MuonResidualsFromTrack:: cov inversion failed!" << std::endl;
0843   return result;
0844 }
0845 
0846 TMatrixD MuonResidualsFromTrack::choleskyCorrMatrix(DetId chamberId) {
0847   bool m_debug = false;
0848 
0849   TMatrixD result(4, 4);
0850   TMatrixDSym corr44 = corrMatrix(chamberId);
0851 
0852   // get an upper triangular matrix U such that corr = U^T * U
0853   TDecompChol decomp(corr44);
0854   bool ok = decomp.Decompose();
0855   result = decomp.GetU();
0856 
0857   if (m_debug)
0858     std::cout << "MuonResidualsFromTrack:: corr cholesky after:" << std::endl;
0859   result.Print();
0860 
0861   if (!ok && m_debug)
0862     std::cout << "MuonResidualsFromTrack:: corr decomposition failed!" << std::endl;
0863   return result;
0864 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!