File indexing completed on 2024-09-07 04:37:47
0001 #include "RecoMuon/MuonIdentification/interface/MuonMesh.h"
0002 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0003 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0004 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0005
0006 #include <utility>
0007 #include <algorithm>
0008
0009 MuonMesh::MuonMesh(const edm::ParameterSet& parm)
0010 : doME1a(parm.getParameter<bool>("ME1a")),
0011 doOverlaps(parm.getParameter<bool>("Overlap")),
0012 doClustering(parm.getParameter<bool>("Clustering")),
0013 OverlapDPhi(parm.getParameter<double>("OverlapDPhi")),
0014 OverlapDTheta(parm.getParameter<double>("OverlapDTheta")),
0015 ClusterDPhi(parm.getParameter<double>("ClusterDPhi")),
0016 ClusterDTheta(parm.getParameter<double>("ClusterDTheta")) {}
0017
0018 void MuonMesh::fillMesh(std::vector<reco::Muon>* inputMuons) {
0019 for (std::vector<reco::Muon>::iterator muonIter1 = inputMuons->begin(); muonIter1 != inputMuons->end(); ++muonIter1) {
0020 if (muonIter1->isTrackerMuon()) {
0021 mesh_[&*muonIter1];
0022 for (std::vector<reco::Muon>::iterator muonIter2 = inputMuons->begin(); muonIter2 != inputMuons->end();
0023 ++muonIter2) {
0024 if (muonIter2->isTrackerMuon()) {
0025 if (muonIter2 != muonIter1) {
0026
0027 for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = muonIter1->matches().begin();
0028 chamberIter1 != muonIter1->matches().end();
0029 ++chamberIter1) {
0030 for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0031 segmentIter1 != chamberIter1->segmentMatches.end();
0032 ++segmentIter1) {
0033 for (std::vector<reco::MuonChamberMatch>::iterator chamberIter2 = muonIter2->matches().begin();
0034 chamberIter2 != muonIter2->matches().end();
0035 ++chamberIter2) {
0036 for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 =
0037 chamberIter2->segmentMatches.begin();
0038 segmentIter2 != chamberIter2->segmentMatches.end();
0039 ++segmentIter2) {
0040
0041
0042 bool addsegment(false);
0043
0044 if (segmentIter1->cscSegmentRef.isNonnull() && segmentIter2->cscSegmentRef.isNonnull()) {
0045 if (doME1a && isDuplicateOf(segmentIter1->cscSegmentRef, segmentIter2->cscSegmentRef) &&
0046 CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(chamberIter2->id).ring() == 4 &&
0047 chamberIter1->id == chamberIter2->id) {
0048 addsegment = true;
0049
0050 }
0051
0052 if (doOverlaps &&
0053 isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0054 std::make_pair(CSCDetId(chamberIter2->id), segmentIter2->cscSegmentRef))) {
0055 addsegment = true;
0056
0057 }
0058
0059 if (doClustering &&
0060 isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0061 std::make_pair(CSCDetId(chamberIter2->id), segmentIter2->cscSegmentRef))) {
0062 addsegment = true;
0063
0064 }
0065
0066 }
0067
0068 if (addsegment) {
0069
0070 if (find(mesh_[&*muonIter1].begin(),
0071 mesh_[&*muonIter1].end(),
0072 std::make_pair(&*muonIter2, std::make_pair(&*chamberIter2, &*segmentIter2))) ==
0073 mesh_[&*muonIter1].end()) {
0074 mesh_[&*muonIter1].push_back(
0075 std::make_pair(&*muonIter2, std::make_pair(&*chamberIter2, &*segmentIter2)));
0076 }
0077 }
0078
0079 }
0080 }
0081 }
0082 }
0083
0084 }
0085 }
0086 }
0087
0088 }
0089 }
0090
0091
0092
0093
0094 if (mesh_.size() == 1) {
0095 for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = mesh_.begin()->first->matches().begin();
0096 chamberIter1 != mesh_.begin()->first->matches().end();
0097 ++chamberIter1) {
0098 for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0099 segmentIter1 != chamberIter1->segmentMatches.end();
0100 ++segmentIter1) {
0101 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0102 }
0103 }
0104 }
0105
0106
0107
0108 if (mesh_.size() > 1) {
0109 for (MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i) {
0110 for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
0111 chamberIter1 != i->first->matches().end();
0112 ++chamberIter1) {
0113 for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0114 segmentIter1 != chamberIter1->segmentMatches.end();
0115 ++segmentIter1) {
0116 bool shared(false);
0117
0118 for (AssociationType::iterator j = i->second.begin(); j != i->second.end(); ++j) {
0119 if (segmentIter1->cscSegmentRef.isNonnull() && j->second.second->cscSegmentRef.isNonnull()) {
0120 if (chamberIter1->id.subdetId() == MuonSubdetId::CSC &&
0121 j->second.first->id.subdetId() == MuonSubdetId::CSC) {
0122 CSCDetId segIterId(chamberIter1->id), shareId(j->second.first->id);
0123
0124 if (doOverlaps && isDuplicateOf(std::make_pair(segIterId, segmentIter1->cscSegmentRef),
0125 std::make_pair(shareId, j->second.second->cscSegmentRef)))
0126 shared = true;
0127
0128 if (doME1a && isDuplicateOf(segmentIter1->cscSegmentRef, j->second.second->cscSegmentRef) &&
0129 segIterId.ring() == 4 && shareId.ring() == 4 && segIterId == segIterId)
0130 shared = true;
0131
0132 if (doClustering &&
0133 isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0134 std::make_pair(CSCDetId(j->second.first->id), j->second.second->cscSegmentRef)))
0135 shared = true;
0136 }
0137 }
0138 }
0139
0140
0141 if (((segmentIter1->mask & 0x1e0000) == 0x1e0000 && !shared) ||
0142 (chamberIter1->id.subdetId() == MuonSubdetId::DT && (segmentIter1->mask & 0x1e0000)))
0143 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0144
0145 }
0146 }
0147 }
0148 }
0149 }
0150
0151 void MuonMesh::pruneMesh() {
0152 for (MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i) {
0153 for (AssociationType::iterator j = i->second.begin(); j != i->second.end(); ++j) {
0154 for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
0155 chamberIter1 != i->first->matches().end();
0156 ++chamberIter1) {
0157 for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0158 segmentIter1 != chamberIter1->segmentMatches.end();
0159 ++segmentIter1) {
0160 if (j->second.second->cscSegmentRef.isNonnull() && segmentIter1->cscSegmentRef.isNonnull()) {
0161
0162
0163
0164 if (doOverlaps &&
0165 isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0166 std::make_pair(CSCDetId(j->second.first->id), j->second.second->cscSegmentRef))) {
0167 if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
0168 j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0169 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
0170 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0171
0172
0173 } else if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ==
0174 j->first->numberOfMatches(
0175 (reco::Muon::ArbitrationType)0x1e0000)) {
0176
0177 if ((segmentIter1->mask & 0x1e0000) >
0178 (j->second.second->mask & 0x1e0000)) {
0179
0180 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
0181 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0182
0183
0184 } else {
0185
0186 }
0187 }
0188 }
0189
0190
0191
0192
0193 if (doME1a && isDuplicateOf(segmentIter1->cscSegmentRef, j->second.second->cscSegmentRef) &&
0194 CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(j->second.first->id).ring() == 4 &&
0195 chamberIter1->id == j->second.first->id) {
0196 if (j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) <
0197 i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0198 for (AssociationType::iterator AsscIter1 = i->second.begin(); AsscIter1 != i->second.end();
0199 ++AsscIter1) {
0200 if (AsscIter1->second.second->cscSegmentRef.isNonnull())
0201 if (j->first == AsscIter1->first && j->second.first == AsscIter1->second.first &&
0202 isDuplicateOf(segmentIter1->cscSegmentRef, AsscIter1->second.second->cscSegmentRef)) {
0203 AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
0204 }
0205 }
0206
0207
0208 } else if (j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ==
0209 i->first->numberOfMatches(
0210 (reco::Muon::ArbitrationType)0x1e0000)) {
0211
0212 bool bestArb(true);
0213
0214 for (AssociationType::iterator AsscIter1 = i->second.begin(); AsscIter1 != i->second.end();
0215 ++AsscIter1) {
0216 if (AsscIter1->second.second->cscSegmentRef.isNonnull())
0217 if (j->first == AsscIter1->first && j->second.first == AsscIter1->second.first &&
0218 isDuplicateOf(segmentIter1->cscSegmentRef, AsscIter1->second.second->cscSegmentRef) &&
0219 (segmentIter1->mask & 0x1e0000) < (AsscIter1->second.second->mask & 0x1e0000))
0220 bestArb = false;
0221 }
0222
0223 if (bestArb) {
0224 for (AssociationType::iterator AsscIter1 = i->second.begin(); AsscIter1 != i->second.end();
0225 ++AsscIter1) {
0226 if (AsscIter1->second.second->cscSegmentRef.isNonnull())
0227 if (j->first == AsscIter1->first && j->second.first == AsscIter1->second.first &&
0228 isDuplicateOf(segmentIter1->cscSegmentRef, AsscIter1->second.second->cscSegmentRef)) {
0229 AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
0230 }
0231 }
0232 }
0233
0234
0235 }
0236 }
0237
0238 if (doClustering &&
0239 isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0240 std::make_pair(CSCDetId(j->second.first->id), j->second.second->cscSegmentRef))) {
0241 if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
0242 j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0243 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0244 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0245
0246
0247 } else if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) <
0248 j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0249 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0250 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0251
0252
0253 } else {
0254
0255 if ((segmentIter1->mask & 0x1e0000) >
0256 (j->second.second->mask & 0x1e0000)) {
0257
0258 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0259 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0260
0261
0262 } else if ((segmentIter1->mask & 0x1e0000) < (j->second.second->mask & 0x1e0000)) {
0263
0264 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0265 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0266
0267
0268 } else {
0269 }
0270 }
0271
0272 }
0273
0274 }
0275 }
0276 }
0277 }
0278 }
0279
0280
0281
0282 for (MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i) {
0283 for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
0284 chamberIter1 != i->first->matches().end();
0285 ++chamberIter1) {
0286 for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0287 segmentIter1 != chamberIter1->segmentMatches.end();
0288 ++segmentIter1) {
0289
0290 if (!segmentIter1->isMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning) && segmentIter1->isMask(0xe00000))
0291 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0292 }
0293 }
0294 }
0295 }
0296
0297 bool MuonMesh::isDuplicateOf(const CSCSegmentRef& lhs, const CSCSegmentRef& rhs)
0298 const
0299 {
0300 bool result(false);
0301
0302 if (!lhs->isME11a_duplicate())
0303 return result;
0304
0305 std::vector<CSCSegment> lhs_duplicates = lhs->duplicateSegments();
0306
0307 if (fabs(lhs->localPosition().x() - rhs->localPosition().x()) < 1E-3 &&
0308 fabs(lhs->localPosition().y() - rhs->localPosition().y()) < 1E-3 &&
0309 fabs(lhs->localDirection().x() / lhs->localDirection().z() -
0310 rhs->localDirection().x() / rhs->localDirection().z()) < 1E-3 &&
0311 fabs(lhs->localDirection().y() / lhs->localDirection().z() -
0312 rhs->localDirection().y() / rhs->localDirection().z()) < 1E-3 &&
0313 fabs(lhs->localPositionError().xx() - rhs->localPositionError().xx()) < 1E-3 &&
0314 fabs(lhs->localPositionError().yy() - rhs->localPositionError().yy()) < 1E-3 &&
0315 fabs(lhs->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
0316 fabs(lhs->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
0317 result = true;
0318
0319 for (std::vector<CSCSegment>::const_iterator segIter1 = lhs_duplicates.begin(); segIter1 != lhs_duplicates.end();
0320 ++segIter1) {
0321
0322 if (fabs(segIter1->localPosition().x() - rhs->localPosition().x()) < 1E-3 &&
0323 fabs(segIter1->localPosition().y() - rhs->localPosition().y()) < 1E-3 &&
0324 fabs(segIter1->localDirection().x() / segIter1->localDirection().z() -
0325 rhs->localDirection().x() / rhs->localDirection().z()) < 1E-3 &&
0326 fabs(segIter1->localDirection().y() / segIter1->localDirection().z() -
0327 rhs->localDirection().y() / rhs->localDirection().z()) < 1E-3 &&
0328 fabs(segIter1->localPositionError().xx() - rhs->localPositionError().xx()) < 1E-3 &&
0329 fabs(segIter1->localPositionError().yy() - rhs->localPositionError().yy()) < 1E-3 &&
0330 fabs(segIter1->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
0331 fabs(segIter1->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
0332 result = true;
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343 }
0344
0345 return result;
0346 }
0347
0348 bool MuonMesh::isDuplicateOf(const std::pair<CSCDetId, CSCSegmentRef>& rhs,
0349 const std::pair<CSCDetId, CSCSegmentRef>& lhs)
0350 const
0351 {
0352 bool result(false);
0353
0354
0355
0356
0357 if (rhs.first.endcap() == lhs.first.endcap() && rhs.first.station() == lhs.first.station() &&
0358 rhs.first.ring() == lhs.first.ring()) {
0359
0360
0361
0362
0363
0364
0365
0366 unsigned modulus = ((rhs.first.ring() != 1 || rhs.first.station() == 1) ? 36 : 18);
0367 int left_neighbor = (((rhs.first.chamber() - 1 + modulus) % modulus == 0)
0368 ? modulus
0369 : (rhs.first.chamber() - 1 + modulus) % modulus);
0370 int right_neighbor = (((rhs.first.chamber() + 1) % modulus == 0) ? modulus : (rhs.first.chamber() + 1) % modulus);
0371
0372 if (lhs.first.chamber() == left_neighbor ||
0373 lhs.first.chamber() == right_neighbor) {
0374
0375 std::vector<CSCSegment> thesegments;
0376 thesegments.push_back(*(lhs.second));
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397 double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
0398 double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();
0399 double rhs_z = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).z();
0400
0401 for (std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413 double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
0414 double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
0415 double lhs_z = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).z();
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430 double phidiff =
0431 (fabs(rhs_phi - lhs_phi) > 2 * M_PI ? fabs(rhs_phi - lhs_phi) - 2 * M_PI : fabs(rhs_phi - lhs_phi));
0432
0433 if (phidiff < OverlapDPhi && fabs(rhs_theta - lhs_theta) < OverlapDTheta && fabs(rhs_z) < fabs(lhs_z) &&
0434 rhs_z * lhs_z > 0)
0435 result = true;
0436 }
0437 }
0438 }
0439
0440 return result;
0441 }
0442
0443 bool MuonMesh::isClusteredWith(const std::pair<CSCDetId, CSCSegmentRef>& lhs,
0444 const std::pair<CSCDetId, CSCSegmentRef>& rhs) const {
0445 bool result(false);
0446
0447
0448
0449
0450
0451
0452 if (rhs.first.endcap() == lhs.first.endcap() && lhs.first.station() < rhs.first.station()) {
0453 std::vector<CSCSegment> thesegments;
0454 thesegments.push_back(*(lhs.second));
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473
0474 double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
0475 double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();
0476
0477 for (std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489 double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
0490 double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505 double phidiff =
0506 (fabs(rhs_phi - lhs_phi) > 2 * M_PI ? fabs(rhs_phi - lhs_phi) - 2 * M_PI : fabs(rhs_phi - lhs_phi));
0507
0508 if (phidiff < ClusterDPhi && fabs(rhs_theta - lhs_theta) < ClusterDTheta)
0509 result = true;
0510 }
0511 }
0512
0513 return result;
0514 }