File indexing completed on 2021-02-14 14:26:17
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #include "RecoMuon/TrackingTools/interface/MuonSegmentMatcher.h"
0013
0014
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
0033 #include <vector>
0034 #include <iostream>
0035
0036 class MuonSegmentMatcher;
0037
0038 using namespace std;
0039
0040
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
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
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
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
0099 for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0100
0101 DetId idT = (*hit)->geographicalId();
0102 if (!(rechit->geographicalId().rawId() == idT.rawId()))
0103 continue;
0104
0105
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
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
0147 if (!(hiti->geographicalId().rawId() == idT.rawId()))
0148 continue;
0149
0150
0151 LocalPoint segLocal = hiti->localPosition();
0152
0153
0154 if ((fabs(pointLocal.x() - segLocal.x()) < ZCutParameter) &&
0155 (fabs(pointLocal.y() - segLocal.y()) < ZCutParameter))
0156 countAgreeingHits++;
0157 }
0158 }
0159
0160 matchRatioZ = countMuonDTHits == 0 ? 0 : countAgreeingHits / countMuonDTHits;
0161 if (nhitsZ)
0162 if (countAgreeingHits / nhitsZ > matchRatioZ)
0163 matchRatioZ = countAgreeingHits / nhitsZ;
0164 }
0165
0166 if (rechit->hasPhi()) {
0167 double countMuonDTHits = 0;
0168 double countAgreeingHits = 0;
0169
0170
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
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();
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
0200 if (!(hiti->geographicalId().rawId() == idT.rawId()))
0201 continue;
0202
0203
0204 LocalPoint segLocal = hiti->localPosition();
0205
0206
0207
0208 if ((fabs(pointLocal.x() - segLocal.x()) < PhiCutParameter) &&
0209 (fabs(pointLocal.y() - segLocal.y()) < PhiCutParameter))
0210 countAgreeingHits++;
0211 }
0212 }
0213
0214 matchRatioPhi = countMuonDTHits != 0 ? countAgreeingHits / countMuonDTHits : 0;
0215 if (nhitsPhi)
0216 if (countAgreeingHits / nhitsPhi > matchRatioPhi)
0217 matchRatioPhi = countAgreeingHits / nhitsPhi;
0218 }
0219
0220 if (dtTightMatch && nhitsPhi && nhitsZ) {
0221 if ((matchRatioPhi > 0.9) && (matchRatioZ > 0.9)) {
0222
0223 pointerTo4DSegments.push_back(&(*rechit));
0224 }
0225 } else {
0226 if ((matchRatioPhi > 0.9 && nhitsPhi) || (matchRatioZ > 0.9 && nhitsZ)) {
0227
0228 pointerTo4DSegments.push_back(&(*rechit));
0229 }
0230 }
0231
0232 }
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
0278 DetId id = hitC->geographicalId();
0279 CSCDetId cscDetIdHit(id.rawId());
0280
0281 if (segments) {
0282 if (!(myChamber.rawId() == cscDetIdHit.rawId()))
0283 continue;
0284
0285
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
0317
0318 if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0319 (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut)) {
0320 CSCcountAgreeingHits++;
0321
0322 }
0323 }
0324 }
0325
0326 matchRatioCSC = countMuonCSCHits == 0 ? 0 : CSCcountAgreeingHits / countMuonCSCHits;
0327
0328 if ((matchRatioCSC > 0.9) && ((countMuonCSCHits > 1) || !cscTightMatch))
0329 pointerToCSCSegments.push_back(&(*segmentCSC));
0330
0331 }
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
0364 DetId id = hitC->geographicalId();
0365 RPCDetId rpcDetIdHit(id.rawId());
0366
0367 if (rpcDetIdHit != myChamber)
0368 continue;
0369 LocalPoint posLocalMuon = hitC->localPosition();
0370
0371
0372
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
0387