File indexing completed on 2024-04-06 12:27:18
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 double CSCXCut = 0.001;
0246 double CSCYCut = 0.001;
0247 double countMuonCSCHits = 0;
0248
0249 for (CSCSegmentCollection::const_iterator segmentCSC = allSegmentsCSC->begin(); segmentCSC != allSegmentsCSC->end();
0250 segmentCSC++) {
0251 double CSCcountAgreeingHits = 0;
0252
0253 if (!segmentCSC->isValid())
0254 continue;
0255
0256 const vector<CSCRecHit2D>& CSCRechits2D = segmentCSC->specificRecHits();
0257 countMuonCSCHits = 0;
0258 CSCDetId myChamber((*segmentCSC).geographicalId().rawId());
0259
0260 bool segments = false;
0261
0262 for (auto const& hitC : muon.recHits()) {
0263 if (!hitC->isValid())
0264 continue;
0265 if (hitC->geographicalId().det() != DetId::Muon)
0266 continue;
0267 if (hitC->geographicalId().subdetId() != MuonSubdetId::CSC)
0268 continue;
0269 if (!hitC->isValid())
0270 continue;
0271 if (hitC->recHits().size() > 1)
0272 segments = true;
0273
0274
0275 DetId id = hitC->geographicalId();
0276 CSCDetId cscDetIdHit(id.rawId());
0277
0278 if (segments) {
0279 if (!(myChamber.rawId() == cscDetIdHit.rawId()))
0280 continue;
0281
0282
0283 LocalPoint positionLocalCSC = hitC->localPosition();
0284 LocalPoint segLocalCSC = segmentCSC->localPosition();
0285 if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0286 (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut))
0287 pointerToCSCSegments.push_back(&(*segmentCSC));
0288 continue;
0289 }
0290
0291 if (!(cscDetIdHit.ring() == myChamber.ring()))
0292 continue;
0293 if (!(cscDetIdHit.station() == myChamber.station()))
0294 continue;
0295 if (!(cscDetIdHit.endcap() == myChamber.endcap()))
0296 continue;
0297 if (!(cscDetIdHit.chamber() == myChamber.chamber()))
0298 continue;
0299
0300 countMuonCSCHits++;
0301
0302 LocalPoint positionLocalCSC = hitC->localPosition();
0303
0304 for (vector<CSCRecHit2D>::const_iterator hiti = CSCRechits2D.begin(); hiti != CSCRechits2D.end(); hiti++) {
0305 if (!hiti->isValid())
0306 continue;
0307 CSCDetId cscDetId((hiti->geographicalId()).rawId());
0308
0309 if (hitC->geographicalId().rawId() != (hiti->geographicalId()).rawId())
0310 continue;
0311
0312 LocalPoint segLocalCSC = hiti->localPosition();
0313
0314
0315 if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0316 (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut)) {
0317 CSCcountAgreeingHits++;
0318
0319 }
0320 }
0321 }
0322
0323 matchRatioCSC = countMuonCSCHits == 0 ? 0 : CSCcountAgreeingHits / countMuonCSCHits;
0324
0325 if ((matchRatioCSC > 0.9) && ((countMuonCSCHits > 1) || !cscTightMatch))
0326 pointerToCSCSegments.push_back(&(*segmentCSC));
0327
0328 }
0329
0330 return pointerToCSCSegments;
0331 }
0332
0333 vector<const RPCRecHit*> MuonSegmentMatcher::matchRPC(const reco::Track& muon, const edm::Event& event) {
0334 using namespace edm;
0335
0336 edm::Handle<RPCRecHitCollection> rpcRecHits;
0337 event.getByToken(rpcRecHitsToken, rpcRecHits);
0338
0339 vector<const RPCRecHit*> pointerToRPCRecHits;
0340 double RPCCut = 0.001;
0341
0342 for (RPCRecHitCollection::const_iterator hitRPC = rpcRecHits->begin(); hitRPC != rpcRecHits->end(); hitRPC++) {
0343 if (!hitRPC->isValid())
0344 continue;
0345
0346 RPCDetId myChamber((*hitRPC).geographicalId().rawId());
0347 LocalPoint posLocalRPC = hitRPC->localPosition();
0348 bool matched = false;
0349
0350 for (auto const& hitC : muon.recHits()) {
0351 if (!hitC->isValid())
0352 continue;
0353 if (hitC->geographicalId().det() != DetId::Muon)
0354 continue;
0355 if (hitC->geographicalId().subdetId() != MuonSubdetId::RPC)
0356 continue;
0357 if (!hitC->isValid())
0358 continue;
0359
0360
0361 DetId id = hitC->geographicalId();
0362 RPCDetId rpcDetIdHit(id.rawId());
0363
0364 if (rpcDetIdHit != myChamber)
0365 continue;
0366 LocalPoint posLocalMuon = hitC->localPosition();
0367
0368
0369
0370 if ((fabs(posLocalMuon.x() - posLocalRPC.x()) < RPCCut)) {
0371 matched = true;
0372 break;
0373 }
0374 }
0375
0376 if (matched)
0377 pointerToRPCRecHits.push_back(&(*hitRPC));
0378 }
0379
0380 return pointerToRPCRecHits;
0381 }
0382
0383
0384