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