TrackingTools/src/MuonSegmentMatcher.cc

0001 // -*- C++ -*- // // Package:  MuonSegmentMatcher // Class:  MuonSegmentMatcher // /**\class MuonSegmentMatcher MuonSegmentMatcher.cc
0002 // Description: <one line class summary>
0003 // Implementation:
0004 //     <Notes on implementation>
0005 //*/
0006 //
0007 // Original Author:  Alan Tua
0008 //         Created:  Wed Jul  9 21:40:17 CEST 2008
0009 //
0010 //
0011
0012 #include "RecoMuon/TrackingTools/interface/MuonSegmentMatcher.h"
0013
0014 // user include files
0015 #include "FWCore/Framework/interface/Event.h"
0016
0017 #include "DataFormats/MuonReco/interface/Muon.h"
0018
0019 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021
0022 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0023 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0024
0025 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0026 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0027 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0028 #include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h"
0029 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0030
0031 // system include files
0032 #include <vector>
0033 #include <iostream>
0034
0035 class MuonSegmentMatcher;
0036
0037 using namespace std;
0038
0039 // constructors and destructor
0040
0041 MuonSegmentMatcher::MuonSegmentMatcher(const edm::ParameterSet& matchParameters, edm::ConsumesCollector& iC)
0042     : DTSegmentTags_(matchParameters.getParameter<edm::InputTag>("DTsegments")),
0043       CSCSegmentTags_(matchParameters.getParameter<edm::InputTag>("CSCsegments")),
0044       dtRadius_(matchParameters.getParameter<double>("DTradius")),
0045       dtTightMatch(matchParameters.getParameter<bool>("TightMatchDT")),
0046       cscTightMatch(matchParameters.getParameter<bool>("TightMatchCSC")) {
0047   if (matchParameters.existsAs<edm::InputTag>("RPChits")) {
0048     RPCHitTags_ = matchParameters.getParameter<edm::InputTag>("RPChits");
0049   } else
0050     RPCHitTags_ = edm::InputTag("hltRpcRecHits");
0051
0052   dtRecHitsToken = iC.consumes<DTRecSegment4DCollection>(DTSegmentTags_);
0053   allSegmentsCSCToken = iC.consumes<CSCSegmentCollection>(CSCSegmentTags_);
0054   rpcRecHitsToken = iC.consumes<RPCRecHitCollection>(RPCHitTags_);
0055 }
0056
0057 void MuonSegmentMatcher::fillPSetDescription(edm::ParameterSetDescription& desc) {
0058   desc.add<edm::InputTag>("DTsegments", edm::InputTag("dt4DSegments"));
0059   desc.add<edm::InputTag>("CSCsegments", edm::InputTag("cscSegments"));
0060   desc.add<double>("DTradius", 0.01);
0061   desc.add<bool>("TightMatchDT", false);
0062   desc.add<bool>("TightMatchCSC", false);
0063   desc.add<edm::InputTag>("RPChits", edm::InputTag("rpcRecHits"));
0064 }
0065
0066 MuonSegmentMatcher::~MuonSegmentMatcher() {}
0067
0068 // member functions
0069 vector<const DTRecSegment4D*> MuonSegmentMatcher::matchDT(const reco::Track& muon, const edm::Event& event) {
0070   using namespace edm;
0071
0072   edm::Handle<DTRecSegment4DCollection> dtRecHits;
0073   event.getByToken(dtRecHitsToken, dtRecHits);
0074
0075   vector<const DTRecSegment4D*> pointerTo4DSegments;
0076
0077   std::vector<TrackingRecHit const*> dtHits;
0078
0079   bool segments = false;
0080
0081   // Loop and select DT recHits
0082   for (auto const& hit : muon.recHits()) {
0083     if (!hit->isValid())
0084       continue;
0085     if (hit->geographicalId().det() != DetId::Muon)
0086       continue;
0087     if (hit->geographicalId().subdetId() != MuonSubdetId::DT)
0088       continue;
0089     if (!hit->recHits().empty())
0090       if ((*hit->recHits().begin())->recHits().size() > 1)
0091         segments = true;
0092     dtHits.push_back(hit);
0093   }
0094
0095   //  cout << "Muon DT hits found: " << dtHits.size() << " segments " << segments << endl;
0096
0097   double PhiCutParameter = dtRadius_;
0098   double ZCutParameter = dtRadius_;
0099   double matchRatioZ = 0;
0100   double matchRatioPhi = 0;
0101
0102   for (DTRecSegment4DCollection::const_iterator rechit = dtRecHits->begin(); rechit != dtRecHits->end(); ++rechit) {
0103     LocalPoint pointLocal = rechit->localPosition();
0104
0105     if (segments) {
0106       // Loop over muon recHits
0107       for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0108         // Pick the one in the same DT Chamber as the muon
0109         DetId idT = (*hit)->geographicalId();
0110         if (!(rechit->geographicalId().rawId() == idT.rawId()))
0111           continue;
0112
0113         // and compare the local positions
0114         LocalPoint segLocal = (*hit)->localPosition();
0115         if ((fabs(pointLocal.x() - segLocal.x()) < ZCutParameter) &&
0116             (fabs(pointLocal.y() - segLocal.y()) < ZCutParameter))
0117           pointerTo4DSegments.push_back(&(*rechit));
0118       }
0119       continue;
0120     }
0121
0122     double nhitsPhi = 0;
0123     double nhitsZ = 0;
0124
0125     if (rechit->hasZed()) {
0126       double countMuonDTHits = 0;
0127       double countAgreeingHits = 0;
0128
0129       const DTRecSegment2D* segmZ;
0130       segmZ = dynamic_cast<const DTRecSegment2D*>(rechit->zSegment());
0131       nhitsZ = segmZ->recHits().size();
0132
0133       const vector<DTRecHit1D> hits1d = segmZ->specificRecHits();
0134       DTChamberId chamberSegIdT((segmZ->geographicalId()).rawId());
0135
0136       // Loop over muon recHits
0137       for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0138         if (!(*hit)->isValid())
0139           continue;
0140
0141         DetId idT = (*hit)->geographicalId();
0142         DTChamberId dtDetIdHitT(idT.rawId());
0143         DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
0144
0145         LocalPoint pointLocal = (*hit)->localPosition();
0146
0147         if ((chamberSegIdT == dtDetIdHitT) && (dtDetLayerIdHitT.superlayer() == 2))
0148           countMuonDTHits++;
0149
0150         for (vector<DTRecHit1D>::const_iterator hiti = hits1d.begin(); hiti != hits1d.end(); hiti++) {
0151           if (!hiti->isValid())
0152             continue;
0153
0154           // Pick the one in the same DT Layer as the 1D hit
0155           if (!(hiti->geographicalId().rawId() == idT.rawId()))
0156             continue;
0157
0158           // and compare the local positions
0159           LocalPoint segLocal = hiti->localPosition();
0160           //      cout << "Zed Segment Point = "<<pointLocal<<"    Muon Point = "<<segLocal<<"  Dist:  "
0161           //           << (fabs(pointLocal.x()-segLocal.x()))+(fabs(pointLocal.y()-segLocal.y()))<< endl;
0162           if ((fabs(pointLocal.x() - segLocal.x()) < ZCutParameter) &&
0163               (fabs(pointLocal.y() - segLocal.y()) < ZCutParameter))
0164             countAgreeingHits++;
0165         }  //End Segment Hit Iteration
0166       }  //End Muon Hit Iteration
0167
0168       matchRatioZ = countMuonDTHits == 0 ? 0 : countAgreeingHits / countMuonDTHits;
0169       if (nhitsZ)
0170         if (countAgreeingHits / nhitsZ > matchRatioZ)
0171           matchRatioZ = countAgreeingHits / nhitsZ;
0172     }  //End HasZed Check
0173
0174     if (rechit->hasPhi()) {
0175       double countMuonDTHits = 0;
0176       double countAgreeingHits = 0;
0177
0178       //PREPARE PARAMETERS FOR SEGMENT DETECTOR GEOMETRY
0179       const DTRecSegment2D* segmPhi;
0180       segmPhi = dynamic_cast<const DTRecSegment2D*>(rechit->phiSegment());
0181       nhitsPhi = segmPhi->recHits().size();
0182
0183       const vector<DTRecHit1D> hits1d = segmPhi->specificRecHits();
0184       DTChamberId chamberSegIdT((segmPhi->geographicalId()).rawId());
0185
0186       // Loop over muon recHits
0187       for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0188         if (!(*hit)->isValid())
0189           continue;
0190
0191         DetId idT = (*hit)->geographicalId();
0192         DTChamberId dtDetIdHitT(idT.rawId());
0193         DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
0194
0195         LocalPoint pointLocal =
0196             (*hit)
0197                 ->localPosition();  //Localposition is in DTLayer http://cmslxr.fnal.gov/lxr/source/DataFormats/DTRecHit/interface/DTRecHit1D.h
0198
0199         if ((chamberSegIdT == dtDetIdHitT) &&
0200             ((dtDetLayerIdHitT.superlayer() == 1) || (dtDetLayerIdHitT.superlayer() == 3)))
0201           countMuonDTHits++;
0202
0203         for (vector<DTRecHit1D>::const_iterator hiti = hits1d.begin(); hiti != hits1d.end(); hiti++) {
0204           if (!hiti->isValid())
0205             continue;
0206
0207           // Pick the one in the same DT Layer as the 1D hit
0208           if (!(hiti->geographicalId().rawId() == idT.rawId()))
0209             continue;
0210
0211           // and compare the local positions
0212           LocalPoint segLocal = hiti->localPosition();
0213           //      cout << "     Phi Segment Point = "<<pointLocal<<"    Muon Point = "<<segLocal<<"  Dist:   "
0214           //           << (fabs(pointLocal.x()-segLocal.x()))+(fabs(pointLocal.y()-segLocal.y()))<< endl;
0215
0216           if ((fabs(pointLocal.x() - segLocal.x()) < PhiCutParameter) &&
0217               (fabs(pointLocal.y() - segLocal.y()) < PhiCutParameter))
0218             countAgreeingHits++;
0219         }  // End Segment Hit Iteration
0220       }  // End Muon Hit Iteration
0221
0222       matchRatioPhi = countMuonDTHits != 0 ? countAgreeingHits / countMuonDTHits : 0;
0223       if (nhitsPhi)
0224         if (countAgreeingHits / nhitsPhi > matchRatioPhi)
0225           matchRatioPhi = countAgreeingHits / nhitsPhi;
0226     }  // End HasPhi Check
0227     //    DTChamberId chamberSegId2((rechit->geographicalId()).rawId());
0228     if (dtTightMatch && nhitsPhi && nhitsZ) {
0229       if ((matchRatioPhi > 0.9) && (matchRatioZ > 0.9)) {
0230         //  cout<<"Making a tight match in Chamber "<<chamberSegId2<<endl;
0231         pointerTo4DSegments.push_back(&(*rechit));
0232       }
0233     } else {
0234       if ((matchRatioPhi > 0.9 && nhitsPhi) || (matchRatioZ > 0.9 && nhitsZ)) {
0235         //  cout<<"Making a loose match in Chamber "<<chamberSegId2<<endl;
0236         pointerTo4DSegments.push_back(&(*rechit));
0237       }
0238     }
0239
0240   }  //End DT Segment Iteration
0241
0242   return pointerTo4DSegments;
0243 }
0244
0245 vector<const CSCSegment*> MuonSegmentMatcher::matchCSC(const reco::Track& muon, const edm::Event& event) {
0246   using namespace edm;
0247
0248   edm::Handle<CSCSegmentCollection> allSegmentsCSC;
0249   event.getByToken(allSegmentsCSCToken, allSegmentsCSC);
0250
0251   vector<const CSCSegment*> pointerToCSCSegments;
0252
0253   double matchRatioCSC = 0;
0254   double CSCXCut = 0.001;
0255   double CSCYCut = 0.001;
0256   double countMuonCSCHits = 0;
0257
0258   for (CSCSegmentCollection::const_iterator segmentCSC = allSegmentsCSC->begin(); segmentCSC != allSegmentsCSC->end();
0259        segmentCSC++) {
0260     double CSCcountAgreeingHits = 0;
0261
0262     if (!segmentCSC->isValid())
0263       continue;
0264
0265     const vector<CSCRecHit2D>& CSCRechits2D = segmentCSC->specificRecHits();
0266     countMuonCSCHits = 0;
0267     CSCDetId myChamber((*segmentCSC).geographicalId().rawId());
0268
0269     bool segments = false;
0270
0271     for (auto const& hitC : muon.recHits()) {
0272       if (!hitC->isValid())
0273         continue;
0274       if (hitC->geographicalId().det() != DetId::Muon)
0275         continue;
0276       if (hitC->geographicalId().subdetId() != MuonSubdetId::CSC)
0277         continue;
0278       if (!hitC->isValid())
0279         continue;
0280       if (hitC->recHits().size() > 1)
0281         segments = true;
0282
0283       //DETECTOR CONSTRUCTION
0284       DetId id = hitC->geographicalId();
0285       CSCDetId cscDetIdHit(id.rawId());
0286
0287       if (segments) {
0288         if (!(myChamber.rawId() == cscDetIdHit.rawId()))
0289           continue;
0290
0291         // and compare the local positions
0292         LocalPoint positionLocalCSC = hitC->localPosition();
0293         LocalPoint segLocalCSC = segmentCSC->localPosition();
0294         if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0295             (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut))
0296           pointerToCSCSegments.push_back(&(*segmentCSC));
0297         continue;
0298       }
0299
0300       if (!(cscDetIdHit.ring() == myChamber.ring()))
0301         continue;
0302       if (!(cscDetIdHit.station() == myChamber.station()))
0303         continue;
0304       if (!(cscDetIdHit.endcap() == myChamber.endcap()))
0305         continue;
0306       if (!(cscDetIdHit.chamber() == myChamber.chamber()))
0307         continue;
0308
0309       countMuonCSCHits++;
0310
0311       LocalPoint positionLocalCSC = hitC->localPosition();
0312
0313       for (vector<CSCRecHit2D>::const_iterator hiti = CSCRechits2D.begin(); hiti != CSCRechits2D.end(); hiti++) {
0314         if (!hiti->isValid())
0315           continue;
0316         CSCDetId cscDetId((hiti->geographicalId()).rawId());
0317
0318         if (hitC->geographicalId().rawId() != (hiti->geographicalId()).rawId())
0319           continue;
0320
0321         LocalPoint segLocalCSC = hiti->localPosition();
0322         //      cout<<"Layer Id (MuonHit) =  "<<cscDetIdHit<<" Muon Local Position (det frame) "<<positionLocalCSC <<endl;
0323         //      cout<<"Layer Id  (CSCHit) =  "<<cscDetId<<"  Hit Local Position (det frame) "<<segLocalCSC <<endl;
0324         if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0325             (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut)) {
0326           CSCcountAgreeingHits++;
0327           //          cout << "   Matched." << endl;
0328         }
0329       }  //End 2D rechit iteration
0330     }  //End muon hit iteration
0331
0332     matchRatioCSC = countMuonCSCHits == 0 ? 0 : CSCcountAgreeingHits / countMuonCSCHits;
0333
0334     if ((matchRatioCSC > 0.9) && ((countMuonCSCHits > 1) || !cscTightMatch))
0335       pointerToCSCSegments.push_back(&(*segmentCSC));
0336
0337   }  //End CSC Segment Iteration
0338
0339   return pointerToCSCSegments;
0340 }
0341
0342 vector<const RPCRecHit*> MuonSegmentMatcher::matchRPC(const reco::Track& muon, const edm::Event& event) {
0343   using namespace edm;
0344
0345   edm::Handle<RPCRecHitCollection> rpcRecHits;
0346   event.getByToken(rpcRecHitsToken, rpcRecHits);
0347
0348   vector<const RPCRecHit*> pointerToRPCRecHits;
0349   double RPCCut = 0.001;
0350
0351   for (RPCRecHitCollection::const_iterator hitRPC = rpcRecHits->begin(); hitRPC != rpcRecHits->end(); hitRPC++) {
0352     if (!hitRPC->isValid())
0353       continue;
0354
0355     RPCDetId myChamber((*hitRPC).geographicalId().rawId());
0356     LocalPoint posLocalRPC = hitRPC->localPosition();
0357     bool matched = false;
0358
0359     for (auto const& hitC : muon.recHits()) {
0360       if (!hitC->isValid())
0361         continue;
0362       if (hitC->geographicalId().det() != DetId::Muon)
0363         continue;
0364       if (hitC->geographicalId().subdetId() != MuonSubdetId::RPC)
0365         continue;
0366       if (!hitC->isValid())
0367         continue;
0368
0369       //DETECTOR CONSTRUCTION
0370       DetId id = hitC->geographicalId();
0371       RPCDetId rpcDetIdHit(id.rawId());
0372
0373       if (rpcDetIdHit != myChamber)
0374         continue;
0375       LocalPoint posLocalMuon = hitC->localPosition();
0376
0377       //        cout<<"Layer Id (MuonHit) =  "<<rpcDetIdHit<<" Muon Local Position (det frame) "<<posLocalMuon <<endl;
0378       //        cout<<"Layer Id  (RPCHit) =  "<<myChamber<<"  Hit Local Position (det frame) "<<posLocalRPC <<endl;
0379       if ((fabs(posLocalMuon.x() - posLocalRPC.x()) < RPCCut)) {
0380         matched = true;
0381         break;
0382       }
0383     }
0384
0385     if (matched)
0386       pointerToRPCRecHits.push_back(&(*hitRPC));
0387   }
0388
0389   return pointerToRPCRecHits;
0390 }
0391
0392 //define this as a plug-in
0393 //DEFINE_FWK_MODULE(MuonSegmentMatcher);