File indexing completed on 2022-10-30 23:36:34
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/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
0032 #include <vector>
0033 #include <iostream>
0034
0035 class MuonSegmentMatcher;
0036
0037 using namespace std;
0038
0039
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 MuonSegmentMatcher::~MuonSegmentMatcher() {}
0058
0059
0060 vector<const DTRecSegment4D*> MuonSegmentMatcher::matchDT(const reco::Track& muon, const edm::Event& event) {
0061 using namespace edm;
0062
0063 edm::Handle<DTRecSegment4DCollection> dtRecHits;
0064 event.getByToken(dtRecHitsToken, dtRecHits);
0065
0066 vector<const DTRecSegment4D*> pointerTo4DSegments;
0067
0068 std::vector<TrackingRecHit const*> dtHits;
0069
0070 bool segments = false;
0071
0072
0073 for (auto const& hit : muon.recHits()) {
0074 if (!hit->isValid())
0075 continue;
0076 if (hit->geographicalId().det() != DetId::Muon)
0077 continue;
0078 if (hit->geographicalId().subdetId() != MuonSubdetId::DT)
0079 continue;
0080 if (!hit->recHits().empty())
0081 if ((*hit->recHits().begin())->recHits().size() > 1)
0082 segments = true;
0083 dtHits.push_back(hit);
0084 }
0085
0086
0087
0088 double PhiCutParameter = dtRadius_;
0089 double ZCutParameter = dtRadius_;
0090 double matchRatioZ = 0;
0091 double matchRatioPhi = 0;
0092
0093 for (DTRecSegment4DCollection::const_iterator rechit = dtRecHits->begin(); rechit != dtRecHits->end(); ++rechit) {
0094 LocalPoint pointLocal = rechit->localPosition();
0095
0096 if (segments) {
0097
0098 for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0099
0100 DetId idT = (*hit)->geographicalId();
0101 if (!(rechit->geographicalId().rawId() == idT.rawId()))
0102 continue;
0103
0104
0105 LocalPoint segLocal = (*hit)->localPosition();
0106 if ((fabs(pointLocal.x() - segLocal.x()) < ZCutParameter) &&
0107 (fabs(pointLocal.y() - segLocal.y()) < ZCutParameter))
0108 pointerTo4DSegments.push_back(&(*rechit));
0109 }
0110 continue;
0111 }
0112
0113 double nhitsPhi = 0;
0114 double nhitsZ = 0;
0115
0116 if (rechit->hasZed()) {
0117 double countMuonDTHits = 0;
0118 double countAgreeingHits = 0;
0119
0120 const DTRecSegment2D* segmZ;
0121 segmZ = dynamic_cast<const DTRecSegment2D*>(rechit->zSegment());
0122 nhitsZ = segmZ->recHits().size();
0123
0124 const vector<DTRecHit1D> hits1d = segmZ->specificRecHits();
0125 DTChamberId chamberSegIdT((segmZ->geographicalId()).rawId());
0126
0127
0128 for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0129 if (!(*hit)->isValid())
0130 continue;
0131
0132 DetId idT = (*hit)->geographicalId();
0133 DTChamberId dtDetIdHitT(idT.rawId());
0134 DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
0135
0136 LocalPoint pointLocal = (*hit)->localPosition();
0137
0138 if ((chamberSegIdT == dtDetIdHitT) && (dtDetLayerIdHitT.superlayer() == 2))
0139 countMuonDTHits++;
0140
0141 for (vector<DTRecHit1D>::const_iterator hiti = hits1d.begin(); hiti != hits1d.end(); hiti++) {
0142 if (!hiti->isValid())
0143 continue;
0144
0145
0146 if (!(hiti->geographicalId().rawId() == idT.rawId()))
0147 continue;
0148
0149
0150 LocalPoint segLocal = hiti->localPosition();
0151
0152
0153 if ((fabs(pointLocal.x() - segLocal.x()) < ZCutParameter) &&
0154 (fabs(pointLocal.y() - segLocal.y()) < ZCutParameter))
0155 countAgreeingHits++;
0156 }
0157 }
0158
0159 matchRatioZ = countMuonDTHits == 0 ? 0 : countAgreeingHits / countMuonDTHits;
0160 if (nhitsZ)
0161 if (countAgreeingHits / nhitsZ > matchRatioZ)
0162 matchRatioZ = countAgreeingHits / nhitsZ;
0163 }
0164
0165 if (rechit->hasPhi()) {
0166 double countMuonDTHits = 0;
0167 double countAgreeingHits = 0;
0168
0169
0170 const DTRecSegment2D* segmPhi;
0171 segmPhi = dynamic_cast<const DTRecSegment2D*>(rechit->phiSegment());
0172 nhitsPhi = segmPhi->recHits().size();
0173
0174 const vector<DTRecHit1D> hits1d = segmPhi->specificRecHits();
0175 DTChamberId chamberSegIdT((segmPhi->geographicalId()).rawId());
0176
0177
0178 for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0179 if (!(*hit)->isValid())
0180 continue;
0181
0182 DetId idT = (*hit)->geographicalId();
0183 DTChamberId dtDetIdHitT(idT.rawId());
0184 DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
0185
0186 LocalPoint pointLocal =
0187 (*hit)
0188 ->localPosition();
0189
0190 if ((chamberSegIdT == dtDetIdHitT) &&
0191 ((dtDetLayerIdHitT.superlayer() == 1) || (dtDetLayerIdHitT.superlayer() == 3)))
0192 countMuonDTHits++;
0193
0194 for (vector<DTRecHit1D>::const_iterator hiti = hits1d.begin(); hiti != hits1d.end(); hiti++) {
0195 if (!hiti->isValid())
0196 continue;
0197
0198
0199 if (!(hiti->geographicalId().rawId() == idT.rawId()))
0200 continue;
0201
0202
0203 LocalPoint segLocal = hiti->localPosition();
0204
0205
0206
0207 if ((fabs(pointLocal.x() - segLocal.x()) < PhiCutParameter) &&
0208 (fabs(pointLocal.y() - segLocal.y()) < PhiCutParameter))
0209 countAgreeingHits++;
0210 }
0211 }
0212
0213 matchRatioPhi = countMuonDTHits != 0 ? countAgreeingHits / countMuonDTHits : 0;
0214 if (nhitsPhi)
0215 if (countAgreeingHits / nhitsPhi > matchRatioPhi)
0216 matchRatioPhi = countAgreeingHits / nhitsPhi;
0217 }
0218
0219 if (dtTightMatch && nhitsPhi && nhitsZ) {
0220 if ((matchRatioPhi > 0.9) && (matchRatioZ > 0.9)) {
0221
0222 pointerTo4DSegments.push_back(&(*rechit));
0223 }
0224 } else {
0225 if ((matchRatioPhi > 0.9 && nhitsPhi) || (matchRatioZ > 0.9 && nhitsZ)) {
0226
0227 pointerTo4DSegments.push_back(&(*rechit));
0228 }
0229 }
0230
0231 }
0232
0233 return pointerTo4DSegments;
0234 }
0235
0236 vector<const CSCSegment*> MuonSegmentMatcher::matchCSC(const reco::Track& muon, const edm::Event& event) {
0237 using namespace edm;
0238
0239 edm::Handle<CSCSegmentCollection> allSegmentsCSC;
0240 event.getByToken(allSegmentsCSCToken, allSegmentsCSC);
0241
0242 vector<const CSCSegment*> pointerToCSCSegments;
0243
0244 double matchRatioCSC = 0;
0245 int numCSC = 0;
0246 double CSCXCut = 0.001;
0247 double CSCYCut = 0.001;
0248 double countMuonCSCHits = 0;
0249
0250 for (CSCSegmentCollection::const_iterator segmentCSC = allSegmentsCSC->begin(); segmentCSC != allSegmentsCSC->end();
0251 segmentCSC++) {
0252 double CSCcountAgreeingHits = 0;
0253
0254 if (!segmentCSC->isValid())
0255 continue;
0256
0257 numCSC++;
0258 const vector<CSCRecHit2D>& CSCRechits2D = segmentCSC->specificRecHits();
0259 countMuonCSCHits = 0;
0260 CSCDetId myChamber((*segmentCSC).geographicalId().rawId());
0261
0262 bool segments = false;
0263
0264 for (auto const& hitC : muon.recHits()) {
0265 if (!hitC->isValid())
0266 continue;
0267 if (hitC->geographicalId().det() != DetId::Muon)
0268 continue;
0269 if (hitC->geographicalId().subdetId() != MuonSubdetId::CSC)
0270 continue;
0271 if (!hitC->isValid())
0272 continue;
0273 if (hitC->recHits().size() > 1)
0274 segments = true;
0275
0276
0277 DetId id = hitC->geographicalId();
0278 CSCDetId cscDetIdHit(id.rawId());
0279
0280 if (segments) {
0281 if (!(myChamber.rawId() == cscDetIdHit.rawId()))
0282 continue;
0283
0284
0285 LocalPoint positionLocalCSC = hitC->localPosition();
0286 LocalPoint segLocalCSC = segmentCSC->localPosition();
0287 if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0288 (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut))
0289 pointerToCSCSegments.push_back(&(*segmentCSC));
0290 continue;
0291 }
0292
0293 if (!(cscDetIdHit.ring() == myChamber.ring()))
0294 continue;
0295 if (!(cscDetIdHit.station() == myChamber.station()))
0296 continue;
0297 if (!(cscDetIdHit.endcap() == myChamber.endcap()))
0298 continue;
0299 if (!(cscDetIdHit.chamber() == myChamber.chamber()))
0300 continue;
0301
0302 countMuonCSCHits++;
0303
0304 LocalPoint positionLocalCSC = hitC->localPosition();
0305
0306 for (vector<CSCRecHit2D>::const_iterator hiti = CSCRechits2D.begin(); hiti != CSCRechits2D.end(); hiti++) {
0307 if (!hiti->isValid())
0308 continue;
0309 CSCDetId cscDetId((hiti->geographicalId()).rawId());
0310
0311 if (hitC->geographicalId().rawId() != (hiti->geographicalId()).rawId())
0312 continue;
0313
0314 LocalPoint segLocalCSC = hiti->localPosition();
0315
0316
0317 if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0318 (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut)) {
0319 CSCcountAgreeingHits++;
0320
0321 }
0322 }
0323 }
0324
0325 matchRatioCSC = countMuonCSCHits == 0 ? 0 : CSCcountAgreeingHits / countMuonCSCHits;
0326
0327 if ((matchRatioCSC > 0.9) && ((countMuonCSCHits > 1) || !cscTightMatch))
0328 pointerToCSCSegments.push_back(&(*segmentCSC));
0329
0330 }
0331
0332 return pointerToCSCSegments;
0333 }
0334
0335 vector<const RPCRecHit*> MuonSegmentMatcher::matchRPC(const reco::Track& muon, const edm::Event& event) {
0336 using namespace edm;
0337
0338 edm::Handle<RPCRecHitCollection> rpcRecHits;
0339 event.getByToken(rpcRecHitsToken, rpcRecHits);
0340
0341 vector<const RPCRecHit*> pointerToRPCRecHits;
0342 double RPCCut = 0.001;
0343
0344 for (RPCRecHitCollection::const_iterator hitRPC = rpcRecHits->begin(); hitRPC != rpcRecHits->end(); hitRPC++) {
0345 if (!hitRPC->isValid())
0346 continue;
0347
0348 RPCDetId myChamber((*hitRPC).geographicalId().rawId());
0349 LocalPoint posLocalRPC = hitRPC->localPosition();
0350 bool matched = false;
0351
0352 for (auto const& hitC : muon.recHits()) {
0353 if (!hitC->isValid())
0354 continue;
0355 if (hitC->geographicalId().det() != DetId::Muon)
0356 continue;
0357 if (hitC->geographicalId().subdetId() != MuonSubdetId::RPC)
0358 continue;
0359 if (!hitC->isValid())
0360 continue;
0361
0362
0363 DetId id = hitC->geographicalId();
0364 RPCDetId rpcDetIdHit(id.rawId());
0365
0366 if (rpcDetIdHit != myChamber)
0367 continue;
0368 LocalPoint posLocalMuon = hitC->localPosition();
0369
0370
0371
0372 if ((fabs(posLocalMuon.x() - posLocalRPC.x()) < RPCCut)) {
0373 matched = true;
0374 break;
0375 }
0376 }
0377
0378 if (matched)
0379 pointerToRPCRecHits.push_back(&(*hitRPC));
0380 }
0381
0382 return pointerToRPCRecHits;
0383 }
0384
0385
0386