Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:32:31 

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