File indexing completed on 2024-09-07 04:35:56
0001 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0004 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0005 #include "DataFormats/VertexReco/interface/Vertex.h"
0006 #include "DataFormats/MuonReco/interface/MuonRPCHitMatch.h"
0007
0008 namespace muon {
0009 SelectionType selectionTypeFromString(const std::string& label) {
0010 SelectionType value = (SelectionType)-1;
0011 bool found = false;
0012 for (int i = 0; selectionTypeStringToEnumMap[i].label && (!found); ++i)
0013 if (!strcmp(label.c_str(), selectionTypeStringToEnumMap[i].label)) {
0014 found = true;
0015 value = selectionTypeStringToEnumMap[i].value;
0016 }
0017
0018
0019 if (!found)
0020 throw cms::Exception("MuonSelectorError") << label << " is not a recognized SelectionType";
0021 return value;
0022 }
0023
0024 reco::Muon::Selector selectorFromString(const std::string& label) {
0025 reco::Muon::Selector value = (reco::Muon::Selector)-1;
0026 bool found = false;
0027 for (int i = 0; selectorStringToEnumMap[i].label && (!found); ++i)
0028 if (!strcmp(label.c_str(), selectorStringToEnumMap[i].label)) {
0029 found = true;
0030 value = selectorStringToEnumMap[i].value;
0031 }
0032
0033
0034 if (!found)
0035 throw cms::Exception("MuonSelectorError") << label << " is not a recognized reco::Muon::Selector";
0036 return value;
0037 }
0038 }
0039
0040 unsigned int muon::RequiredStationMask(const reco::Muon& muon,
0041 double maxChamberDist,
0042 double maxChamberDistPull,
0043 reco::Muon::ArbitrationType arbitrationType) {
0044 unsigned int theMask = 0;
0045
0046 for (int stationIdx = 1; stationIdx < 5; ++stationIdx) {
0047 for (int detectorIdx = 1; detectorIdx < 3; ++detectorIdx) {
0048 float dist = muon.trackDist(stationIdx, detectorIdx, arbitrationType);
0049 if (dist < maxChamberDist &&
0050 dist < maxChamberDistPull * muon.trackDistErr(stationIdx, detectorIdx, arbitrationType))
0051 theMask += 1 << ((stationIdx - 1) + 4 * (detectorIdx - 1));
0052 }
0053 }
0054 return theMask;
0055 }
0056
0057
0058 float muon::caloCompatibility(const reco::Muon& muon) { return muon.caloCompatibility(); }
0059
0060
0061 float muon::segmentCompatibility(const reco::Muon& muon, reco::Muon::ArbitrationType arbitrationType) {
0062 bool use_weight_regain_at_chamber_boundary = true;
0063 bool use_match_dist_penalty = true;
0064
0065 int nr_of_stations_crossed = 0;
0066 std::vector<int> stations_w_track(8);
0067 std::vector<int> station_has_segmentmatch(8);
0068 std::vector<int> station_was_crossed(8);
0069 std::vector<float> stations_w_track_at_boundary(8);
0070 std::vector<float> station_weight(8);
0071 int position_in_stations = 0;
0072 float full_weight = 0.;
0073
0074 for (int i = 1; i <= 8; ++i) {
0075
0076
0077
0078 if (i <= 4) {
0079 float thisTrackDist = muon.trackDist(i, 1, arbitrationType);
0080 if (thisTrackDist < 999999) {
0081 ++nr_of_stations_crossed;
0082 station_was_crossed[i - 1] = 1;
0083 if (thisTrackDist > -10.)
0084 stations_w_track_at_boundary[i - 1] = thisTrackDist;
0085 else
0086 stations_w_track_at_boundary[i - 1] = 0.;
0087 }
0088
0089 if (muon.segmentX(i, 1, arbitrationType) < 999999) {
0090 station_has_segmentmatch[i - 1] = 1;
0091 }
0092 } else {
0093 float thisTrackDist = muon.trackDist(i - 4, 2, arbitrationType);
0094 if (thisTrackDist < 999999) {
0095 ++nr_of_stations_crossed;
0096 station_was_crossed[i - 1] = 1;
0097 if (thisTrackDist > -10.)
0098 stations_w_track_at_boundary[i - 1] = thisTrackDist;
0099 else
0100 stations_w_track_at_boundary[i - 1] = 0.;
0101 }
0102
0103 if (muon.segmentX(i - 4, 2, arbitrationType) < 999999) {
0104 station_has_segmentmatch[i - 1] = 1;
0105 }
0106 }
0107
0108
0109
0110
0111
0112
0113
0114
0115 }
0116
0117
0118
0119
0120
0121
0122
0123 const float attenuate_weight_regain = 0.5;
0124
0125 for (int i = 1; i <= 8; ++i) {
0126
0127
0128
0129
0130
0131 if (station_was_crossed[i - 1] > 0) {
0132
0133
0134
0135
0136
0137 ++position_in_stations;
0138
0139 switch (nr_of_stations_crossed) {
0140 case 1:
0141 station_weight[i - 1] = 1.f;
0142 break;
0143 case 2:
0144 if (position_in_stations == 1)
0145 station_weight[i - 1] = 0.33f;
0146 else
0147 station_weight[i - 1] = 0.67f;
0148 break;
0149 case 3:
0150 if (position_in_stations == 1)
0151 station_weight[i - 1] = 0.23f;
0152 else if (position_in_stations == 2)
0153 station_weight[i - 1] = 0.33f;
0154 else
0155 station_weight[i - 1] = 0.44f;
0156 break;
0157 case 4:
0158 if (position_in_stations == 1)
0159 station_weight[i - 1] = 0.10f;
0160 else if (position_in_stations == 2)
0161 station_weight[i - 1] = 0.20f;
0162 else if (position_in_stations == 3)
0163 station_weight[i - 1] = 0.30f;
0164 else
0165 station_weight[i - 1] = 0.40f;
0166 break;
0167
0168 default:
0169
0170
0171
0172 station_weight[i - 1] = 1.f / nr_of_stations_crossed;
0173 }
0174
0175 if (use_weight_regain_at_chamber_boundary) {
0176 if (station_has_segmentmatch[i - 1] <= 0 && stations_w_track_at_boundary[i - 1] != 0.) {
0177
0178
0179
0180
0181 station_weight[i - 1] = station_weight[i - 1] * attenuate_weight_regain * 0.5f *
0182 (std::erf(stations_w_track_at_boundary[i - 1] / 6.f) + 1.f);
0183 } else if (station_has_segmentmatch[i - 1] <= 0 &&
0184 stations_w_track_at_boundary[i - 1] == 0.f) {
0185
0186 station_weight[i - 1] = 0.f;
0187 }
0188 } else {
0189 if (station_has_segmentmatch[i - 1] <= 0)
0190 station_weight[i - 1] = 0.f;
0191 }
0192
0193
0194 if (station_has_segmentmatch[i - 1] > 0) {
0195 if (i <= 4) {
0196 if (muon.dY(i, 1, arbitrationType) < 999999.f &&
0197 muon.dX(i, 1, arbitrationType) < 999999.f) {
0198 const float pullTot2 =
0199 std::pow(muon.pullX(i, 1, arbitrationType), 2.) + std::pow(muon.pullY(i, 1, arbitrationType), 2.);
0200 if (pullTot2 > 1.f) {
0201 const float dxy2 =
0202 std::pow(muon.dX(i, 1, arbitrationType), 2.) + std::pow(muon.dY(i, 1, arbitrationType), 2.);
0203
0204 if (use_match_dist_penalty) {
0205
0206 if (dxy2 < 9.f && pullTot2 > 9.f) {
0207 if (dxy2 > 1.f)
0208 station_weight[i - 1] *= 1.f / std::pow(dxy2, .125);
0209 } else {
0210 station_weight[i - 1] *= 1.f / std::pow(pullTot2, .125);
0211 }
0212 }
0213 }
0214 } else if (muon.dY(i, 1, arbitrationType) >= 999999.f) {
0215
0216 if (muon.pullX(i, 1, arbitrationType) > 1.f) {
0217
0218 if (use_match_dist_penalty) {
0219
0220 if (muon.dX(i, 1, arbitrationType) < 3.f && muon.pullX(i, 1, arbitrationType) > 3.f) {
0221 if (muon.dX(i, 1, arbitrationType) > 1.f)
0222 station_weight[i - 1] *= 1.f / std::pow(muon.dX(i, 1, arbitrationType), .25);
0223 } else {
0224 station_weight[i - 1] *= 1.f / std::pow(muon.pullX(i, 1, arbitrationType), .25);
0225 }
0226 }
0227 }
0228 } else {
0229
0230 if (muon.pullY(i, 1, arbitrationType) > 1.f) {
0231
0232 if (use_match_dist_penalty) {
0233
0234 if (muon.dY(i, 1, arbitrationType) < 3. && muon.pullY(i, 1, arbitrationType) > 3.) {
0235 if (muon.dY(i, 1, arbitrationType) > 1.f)
0236 station_weight[i - 1] *= 1.f / std::pow(muon.dY(i, 1, arbitrationType), .25);
0237 } else {
0238 station_weight[i - 1] *= 1.f / std::pow(muon.pullY(i, 1, arbitrationType), .25);
0239 }
0240 }
0241 }
0242 }
0243 } else {
0244 const float pullTot2 =
0245 std::pow(muon.pullX(i - 4, 2, arbitrationType), 2.) + std::pow(muon.pullY(i - 4, 2, arbitrationType), 2.);
0246 if (pullTot2 > 1.f) {
0247
0248 if (use_match_dist_penalty) {
0249 const float dxy2 =
0250 std::pow(muon.dX(i - 4, 2, arbitrationType), 2.) + std::pow(muon.dY(i - 4, 2, arbitrationType), 2.);
0251
0252 if (dxy2 < 9.f && pullTot2 > 9.f) {
0253 if (dxy2 > 1.f)
0254 station_weight[i - 1] *= 1.f / std::pow(dxy2, .125);
0255 } else {
0256 station_weight[i - 1] *= 1.f / std::pow(pullTot2, .125);
0257 }
0258 }
0259 }
0260 }
0261 }
0262
0263
0264
0265
0266
0267 } else {
0268 station_weight[i - 1] = 0.f;
0269 }
0270
0271
0272 full_weight += station_weight[i - 1];
0273 }
0274
0275
0276
0277
0278 if (nr_of_stations_crossed == 0) {
0279
0280 full_weight = 0.5f;
0281 }
0282
0283
0284
0285
0286
0287 return full_weight;
0288 }
0289
0290 bool muon::isGoodMuon(const reco::Muon& muon,
0291 AlgorithmType type,
0292 double minCompatibility,
0293 reco::Muon::ArbitrationType arbitrationType) {
0294 if (!muon.isMatchesValid())
0295 return false;
0296 bool goodMuon = false;
0297
0298 switch (type) {
0299 case TM2DCompatibility:
0300
0301 if (((0.8 * caloCompatibility(muon)) + (1.2 * segmentCompatibility(muon, arbitrationType))) > minCompatibility)
0302 goodMuon = true;
0303 else
0304 goodMuon = false;
0305 return goodMuon;
0306 break;
0307 default:
0308
0309 goodMuon = false;
0310 return goodMuon;
0311 break;
0312 }
0313 }
0314
0315 bool muon::isGoodMuon(const reco::Muon& muon,
0316 AlgorithmType type,
0317 int minNumberOfMatches,
0318 double maxAbsDx,
0319 double maxAbsPullX,
0320 double maxAbsDy,
0321 double maxAbsPullY,
0322 double maxChamberDist,
0323 double maxChamberDistPull,
0324 reco::Muon::ArbitrationType arbitrationType,
0325 bool syncMinNMatchesNRequiredStationsInBarrelOnly,
0326 bool applyAlsoAngularCuts) {
0327 if (!muon.isMatchesValid())
0328 return false;
0329 bool goodMuon = false;
0330
0331 if (type == TMLastStation) {
0332
0333
0334 if (minNumberOfMatches == 0)
0335 return true;
0336
0337 unsigned int theStationMask = muon.stationMask(arbitrationType);
0338 unsigned int theRequiredStationMask =
0339 RequiredStationMask(muon, maxChamberDist, maxChamberDistPull, arbitrationType);
0340
0341
0342 int numSegs = 0;
0343 int numRequiredStations = 0;
0344 for (int it = 0; it < 8; ++it) {
0345 if (theStationMask & 1 << it)
0346 ++numSegs;
0347 if (theRequiredStationMask & 1 << it)
0348 ++numRequiredStations;
0349 }
0350
0351
0352
0353 if (syncMinNMatchesNRequiredStationsInBarrelOnly) {
0354
0355 if (std::abs(muon.eta()) < 1.2) {
0356 if (minNumberOfMatches > numRequiredStations)
0357 minNumberOfMatches = numRequiredStations;
0358 if (minNumberOfMatches < 1)
0359 minNumberOfMatches = 1;
0360 }
0361 } else {
0362 if (minNumberOfMatches > numRequiredStations)
0363 minNumberOfMatches = numRequiredStations;
0364 if (minNumberOfMatches < 1)
0365 minNumberOfMatches = 1;
0366 }
0367
0368 if (numSegs >= minNumberOfMatches)
0369 goodMuon = true;
0370
0371
0372
0373
0374
0375 int lastSegBit = 0;
0376 if (theRequiredStationMask) {
0377 for (int stationIdx = 7; stationIdx >= 0; --stationIdx)
0378 if (theRequiredStationMask & 1 << stationIdx) {
0379 if (theStationMask & 1 << stationIdx) {
0380 lastSegBit = stationIdx;
0381 goodMuon &= 1;
0382 break;
0383 } else {
0384 goodMuon = false;
0385 break;
0386 }
0387 }
0388 } else {
0389 for (int stationIdx = 7; stationIdx >= 0; --stationIdx)
0390 if (theStationMask & 1 << stationIdx) {
0391 lastSegBit = stationIdx;
0392 break;
0393 }
0394 }
0395
0396 if (!goodMuon)
0397 return false;
0398
0399
0400 int station = 0, detector = 0;
0401 station = lastSegBit < 4 ? lastSegBit + 1 : lastSegBit - 3;
0402 detector = lastSegBit < 4 ? 1 : 2;
0403
0404
0405 if (std::abs(muon.pullX(station, detector, arbitrationType, true)) > maxAbsPullX &&
0406 std::abs(muon.dX(station, detector, arbitrationType)) > maxAbsDx)
0407 return false;
0408
0409 if (applyAlsoAngularCuts && std::abs(muon.pullDxDz(station, detector, arbitrationType, true)) > maxAbsPullX)
0410 return false;
0411
0412
0413 if (maxAbsDy < 999999) {
0414
0415
0416 if (detector == 2) {
0417 if (std::abs(muon.pullY(station, 2, arbitrationType, true)) > maxAbsPullY &&
0418 std::abs(muon.dY(station, 2, arbitrationType)) > maxAbsDy)
0419 return false;
0420
0421 if (applyAlsoAngularCuts && std::abs(muon.pullDyDz(station, 2, arbitrationType, true)) > maxAbsPullY)
0422 return false;
0423 } else {
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439 for (int stationIdx = station; stationIdx > 0; --stationIdx) {
0440 if (!(theStationMask & 1 << (stationIdx - 1)))
0441 continue;
0442
0443 if (muon.dY(stationIdx, 1, arbitrationType) > 999998)
0444 continue;
0445
0446 if (std::abs(muon.pullY(stationIdx, 1, arbitrationType, true)) > maxAbsPullY &&
0447 std::abs(muon.dY(stationIdx, 1, arbitrationType)) > maxAbsDy) {
0448 return false;
0449 }
0450
0451 if (applyAlsoAngularCuts && std::abs(muon.pullDyDz(stationIdx, 1, arbitrationType, true)) > maxAbsPullY)
0452 return false;
0453
0454
0455 return true;
0456 }
0457 }
0458 }
0459
0460 return goodMuon;
0461 }
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473 if (type == TMOneStation) {
0474 unsigned int theStationMask = muon.stationMask(arbitrationType);
0475
0476
0477 if (!theStationMask)
0478 return false;
0479
0480 int station = 0, detector = 0;
0481
0482
0483
0484
0485 bool existsGoodDTSegX = false;
0486 bool existsDTSegY = false;
0487
0488
0489
0490 for (int stationIdx = 0; stationIdx <= 7; ++stationIdx)
0491 if (theStationMask & 1 << stationIdx) {
0492 station = stationIdx < 4 ? stationIdx + 1 : stationIdx - 3;
0493 detector = stationIdx < 4 ? 1 : 2;
0494
0495 if ((std::abs(muon.pullX(station, detector, arbitrationType, true)) > maxAbsPullX &&
0496 std::abs(muon.dX(station, detector, arbitrationType)) > maxAbsDx) ||
0497 (applyAlsoAngularCuts && std::abs(muon.pullDxDz(station, detector, arbitrationType, true)) > maxAbsPullX))
0498 continue;
0499 else if (detector == 1)
0500 existsGoodDTSegX = true;
0501
0502
0503 if (maxAbsDy < 999999) {
0504 if (detector == 2) {
0505 if ((std::abs(muon.pullY(station, 2, arbitrationType, true)) > maxAbsPullY &&
0506 std::abs(muon.dY(station, 2, arbitrationType)) > maxAbsDy) ||
0507 (applyAlsoAngularCuts && std::abs(muon.pullDyDz(station, 2, arbitrationType, true)) > maxAbsPullY))
0508 continue;
0509 } else {
0510 if (muon.dY(station, 1, arbitrationType) > 999998)
0511 continue;
0512 else
0513 existsDTSegY = true;
0514
0515 if ((std::abs(muon.pullY(station, 1, arbitrationType, true)) > maxAbsPullY &&
0516 std::abs(muon.dY(station, 1, arbitrationType)) > maxAbsDy) ||
0517 (applyAlsoAngularCuts && std::abs(muon.pullDyDz(station, 1, arbitrationType, true)) > maxAbsPullY)) {
0518 continue;
0519 }
0520 }
0521 }
0522
0523
0524 return true;
0525 }
0526
0527
0528
0529
0530
0531
0532
0533
0534 if (maxAbsDy < 999999) {
0535 if (existsDTSegY)
0536 return false;
0537 else if (existsGoodDTSegX)
0538 return true;
0539 } else
0540 return false;
0541 }
0542
0543 if (type == RPCMu) {
0544 if (minNumberOfMatches == 0)
0545 return true;
0546
0547 int nMatch = 0;
0548 for (const auto& chamberMatch : muon.matches()) {
0549 if (chamberMatch.detector() != MuonSubdetId::RPC)
0550 continue;
0551
0552 const float trkX = chamberMatch.x;
0553 const float errX = chamberMatch.xErr;
0554
0555 for (const auto& rpcMatch : chamberMatch.rpcMatches) {
0556 const float rpcX = rpcMatch.x;
0557 const float dX = std::abs(rpcX - trkX);
0558 if (dX < maxAbsDx or dX < maxAbsPullX * errX) {
0559 ++nMatch;
0560 break;
0561 }
0562 }
0563 }
0564
0565 if (nMatch >= minNumberOfMatches)
0566 return true;
0567 else
0568 return false;
0569 }
0570
0571 if (type == ME0Mu) {
0572 if (minNumberOfMatches == 0)
0573 return true;
0574
0575 int nMatch = 0;
0576 for (const auto& chamberMatch : muon.matches()) {
0577 if (chamberMatch.detector() != MuonSubdetId::ME0)
0578 continue;
0579
0580 const float trkX = chamberMatch.x;
0581 const float errX2 = chamberMatch.xErr * chamberMatch.xErr;
0582 const float trkY = chamberMatch.y;
0583 const float errY2 = chamberMatch.yErr * chamberMatch.yErr;
0584
0585 for (const auto& segment : chamberMatch.me0Matches) {
0586 const float me0X = segment.x;
0587 const float me0ErrX2 = segment.xErr * segment.xErr;
0588 const float me0Y = segment.y;
0589 const float me0ErrY2 = segment.yErr * segment.yErr;
0590
0591 const float dX = std::abs(me0X - trkX);
0592 const float dY = std::abs(me0Y - trkY);
0593 const float invPullX2 = errX2 + me0ErrX2;
0594 const float invPullY2 = errY2 + me0ErrY2;
0595
0596 if ((dX < maxAbsDx or dX < maxAbsPullX * std::sqrt(invPullX2)) and
0597 (dY < maxAbsDy or dY < maxAbsPullY * std::sqrt(invPullY2))) {
0598 ++nMatch;
0599 break;
0600 }
0601 }
0602 }
0603
0604 return (nMatch >= minNumberOfMatches);
0605 }
0606
0607 if (type == GEMMu) {
0608 if (minNumberOfMatches == 0)
0609 return true;
0610
0611 int nMatch = 0;
0612 for (const auto& chamberMatch : muon.matches()) {
0613 if (chamberMatch.detector() != MuonSubdetId::GEM)
0614 continue;
0615
0616 const float trkX = chamberMatch.x;
0617 const float errX2 = chamberMatch.xErr * chamberMatch.xErr;
0618 const float trkY = chamberMatch.y;
0619 const float errY2 = chamberMatch.yErr * chamberMatch.yErr;
0620
0621 for (const auto& segment : chamberMatch.gemMatches) {
0622 const float gemX = segment.x;
0623 const float gemErrX2 = segment.xErr * segment.xErr;
0624 const float gemY = segment.y;
0625 const float gemErrY2 = segment.yErr * segment.yErr;
0626
0627 const float dX = std::abs(gemX - trkX);
0628 const float dY = std::abs(gemY - trkY);
0629 const float invPullX2 = errX2 + gemErrX2;
0630 const float invPullY2 = errY2 + gemErrY2;
0631
0632 if ((dX < maxAbsDx or dX < maxAbsPullX * std::sqrt(invPullX2)) and
0633 (dY < maxAbsDy or dY < maxAbsPullY * std::sqrt(invPullY2))) {
0634 ++nMatch;
0635 break;
0636 }
0637 }
0638 }
0639
0640 return (nMatch >= minNumberOfMatches);
0641 }
0642
0643 return goodMuon;
0644 }
0645
0646 bool muon::isGoodMuon(const reco::Muon& muon, SelectionType type, reco::Muon::ArbitrationType arbitrationType) {
0647 switch (type) {
0648 case muon::All:
0649 return true;
0650 break;
0651 case muon::AllGlobalMuons:
0652 return muon.isGlobalMuon();
0653 break;
0654 case muon::AllTrackerMuons:
0655 return muon.isTrackerMuon();
0656 break;
0657 case muon::AllStandAloneMuons:
0658 return muon.isStandAloneMuon();
0659 break;
0660 case muon::TrackerMuonArbitrated:
0661 return muon.isTrackerMuon() && muon.numberOfMatches(arbitrationType) > 0;
0662 break;
0663 case muon::AllArbitrated:
0664 return !muon.isTrackerMuon() || muon.numberOfMatches(arbitrationType) > 0;
0665 break;
0666 case muon::GlobalMuonPromptTight:
0667 return muon.isGlobalMuon() && muon.globalTrack()->normalizedChi2() < 10. &&
0668 muon.globalTrack()->hitPattern().numberOfValidMuonHits() > 0;
0669 break;
0670
0671
0672
0673
0674
0675
0676
0677
0678
0679
0680 case muon::TMLastStationLoose:
0681 return muon.isTrackerMuon() &&
0682 isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, true, false);
0683 break;
0684 case muon::TMLastStationTight:
0685 return muon.isTrackerMuon() &&
0686 isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, true, false);
0687 break;
0688 case muon::TMOneStationLoose:
0689 return muon.isTrackerMuon() &&
0690 isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, false);
0691 break;
0692 case muon::TMOneStationTight:
0693 return muon.isTrackerMuon() &&
0694 isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, false);
0695 break;
0696 case muon::TMLastStationOptimizedLowPtLoose:
0697 if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
0698 return muon.isTrackerMuon() &&
0699 isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, false);
0700 else
0701 return muon.isTrackerMuon() &&
0702 isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, false, false);
0703 break;
0704 case muon::TMLastStationOptimizedLowPtTight:
0705 if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
0706 return muon.isTrackerMuon() &&
0707 isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, false);
0708 else
0709 return muon.isTrackerMuon() &&
0710 isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, false, false);
0711 break;
0712
0713 case muon::TM2DCompatibilityLoose:
0714 return muon.isTrackerMuon() && isGoodMuon(muon, TM2DCompatibility, 0.7, arbitrationType);
0715 break;
0716
0717 case muon::TM2DCompatibilityTight:
0718 return muon.isTrackerMuon() && isGoodMuon(muon, TM2DCompatibility, 1.0, arbitrationType);
0719 break;
0720 case muon::GMTkChiCompatibility:
0721 return muon.isGlobalMuon() && muon.isQualityValid() &&
0722 std::abs(muon.combinedQuality().trkRelChi2 - muon.innerTrack()->normalizedChi2()) < 2.0;
0723 break;
0724 case muon::GMStaChiCompatibility:
0725 return muon.isGlobalMuon() && muon.isQualityValid() &&
0726 std::abs(muon.combinedQuality().staRelChi2 - muon.outerTrack()->normalizedChi2()) < 2.0;
0727 break;
0728 case muon::GMTkKinkTight:
0729 return muon.isGlobalMuon() && muon.isQualityValid() && muon.combinedQuality().trkKink < 100.0;
0730 break;
0731 case muon::TMLastStationAngLoose:
0732 return muon.isTrackerMuon() &&
0733 isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, false, true);
0734 break;
0735 case muon::TMLastStationAngTight:
0736 return muon.isTrackerMuon() &&
0737 isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, false, true);
0738 break;
0739 case muon::TMOneStationAngLoose:
0740 return muon.isTrackerMuon() &&
0741 isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, true);
0742 break;
0743 case muon::TMOneStationAngTight:
0744 return muon.isTrackerMuon() &&
0745 isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, true);
0746 break;
0747 case muon::TMLastStationOptimizedBarrelLowPtLoose:
0748 if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
0749 return muon.isTrackerMuon() &&
0750 isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, false);
0751 else
0752 return muon.isTrackerMuon() &&
0753 isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, true, false);
0754 break;
0755 case muon::TMLastStationOptimizedBarrelLowPtTight:
0756 if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
0757 return muon.isTrackerMuon() &&
0758 isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, false);
0759 else
0760 return muon.isTrackerMuon() &&
0761 isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, true, false);
0762 break;
0763 case muon::RPCMuLoose:
0764 return muon.isRPCMuon() && isGoodMuon(muon, RPCMu, 2, 20, 4, 1e9, 1e9, 1e9, 1e9, arbitrationType, false, false);
0765 break;
0766 case muon::AllME0Muons:
0767 return muon.isME0Muon();
0768 break;
0769 case muon::ME0MuonArbitrated:
0770 return muon.isME0Muon() &&
0771 isGoodMuon(muon, ME0Mu, 1, 1e9, 1e9, 1e9, 1e9, 1e9, 1e9, arbitrationType, false, false);
0772 break;
0773 case muon::AllGEMMuons:
0774 return muon.isGEMMuon();
0775 break;
0776 case muon::GEMMuonArbitrated:
0777 return muon.isGEMMuon() &&
0778 isGoodMuon(muon, GEMMu, 1, 1e9, 1e9, 1e9, 1e9, 1e9, 1e9, arbitrationType, false, false);
0779 break;
0780 case muon::TriggerIdLoose:
0781 return isLooseTriggerMuon(muon);
0782 break;
0783 default:
0784 return false;
0785 }
0786 }
0787
0788 bool muon::overlap(
0789 const reco::Muon& muon1, const reco::Muon& muon2, double pullX, double pullY, bool checkAdjacentChambers) {
0790 unsigned int nMatches1 = muon1.numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
0791 unsigned int nMatches2 = muon2.numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
0792 unsigned int betterMuon = (muon1.pt() > muon2.pt() ? 1 : 2);
0793 for (std::vector<reco::MuonChamberMatch>::const_iterator chamber1 = muon1.matches().begin();
0794 chamber1 != muon1.matches().end();
0795 ++chamber1)
0796 for (std::vector<reco::MuonChamberMatch>::const_iterator chamber2 = muon2.matches().begin();
0797 chamber2 != muon2.matches().end();
0798 ++chamber2) {
0799
0800
0801
0802
0803 if (chamber1->id == chamber2->id) {
0804
0805 if (std::abs(chamber1->x - chamber2->x) <
0806 pullX * sqrt(chamber1->xErr * chamber1->xErr + chamber2->xErr * chamber2->xErr)) {
0807 if (betterMuon == 1)
0808 nMatches2--;
0809 else
0810 nMatches1--;
0811 if (nMatches1 == 0 || nMatches2 == 0)
0812 return true;
0813 continue;
0814 }
0815 if (std::abs(chamber1->y - chamber2->y) <
0816 pullY * sqrt(chamber1->yErr * chamber1->yErr + chamber2->yErr * chamber2->yErr)) {
0817 if (betterMuon == 1)
0818 nMatches2--;
0819 else
0820 nMatches1--;
0821 if (nMatches1 == 0 || nMatches2 == 0)
0822 return true;
0823 }
0824 } else {
0825 if (!checkAdjacentChambers)
0826 continue;
0827
0828 if (chamber1->id.subdetId() != MuonSubdetId::CSC || chamber2->id.subdetId() != MuonSubdetId::CSC)
0829 continue;
0830 CSCDetId id1(chamber1->id);
0831 CSCDetId id2(chamber2->id);
0832 if (id1.endcap() != id2.endcap())
0833 continue;
0834 if (id1.station() != id2.station())
0835 continue;
0836 if (id1.ring() != id2.ring())
0837 continue;
0838 if (std::abs(id1.chamber() - id2.chamber()) > 1)
0839 continue;
0840
0841
0842
0843
0844
0845
0846 if (std::abs(chamber1->edgeX) > chamber1->xErr * pullX)
0847 continue;
0848 if (std::abs(chamber2->edgeX) > chamber2->xErr * pullX)
0849 continue;
0850 if (chamber1->x * chamber2->x < 0) {
0851 if (betterMuon == 1)
0852 nMatches2--;
0853 else
0854 nMatches1--;
0855 if (nMatches1 == 0 || nMatches2 == 0)
0856 return true;
0857 }
0858 }
0859 }
0860 return false;
0861 }
0862
0863 bool muon::isLooseTriggerMuon(const reco::Muon& muon) {
0864
0865
0866
0867 bool tk_id = muon::isGoodMuon(muon, TMOneStationTight);
0868 if (not tk_id)
0869 return false;
0870 bool layer_requirements = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5;
0871 bool match_requirements =
0872 (muon.expectedNnumberOfMatchedStations() < 2) or (muon.numberOfMatchedStations() > 1) or (muon.pt() < 8);
0873 return layer_requirements and match_requirements;
0874 }
0875
0876 bool muon::isTightMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
0877 if (!muon.isPFMuon() || !muon.isGlobalMuon())
0878 return false;
0879
0880 bool muID = isGoodMuon(muon, GlobalMuonPromptTight) && (muon.numberOfMatchedStations() > 1);
0881
0882 bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
0883 muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0;
0884
0885 bool ip = std::abs(muon.muonBestTrack()->dxy(vtx.position())) < 0.2 &&
0886 std::abs(muon.muonBestTrack()->dz(vtx.position())) < 0.5;
0887
0888 return muID && hits && ip;
0889 }
0890
0891 bool muon::isLooseMuon(const reco::Muon& muon) {
0892 return muon.isPFMuon() && (muon.isGlobalMuon() || muon.isTrackerMuon());
0893 }
0894
0895 bool muon::isMediumMuon(const reco::Muon& muon, bool run2016_hip_mitigation) {
0896 if (not isLooseMuon(muon))
0897 return false;
0898 if (run2016_hip_mitigation) {
0899 if (muon.innerTrack()->validFraction() < 0.49)
0900 return false;
0901 } else {
0902 if (muon.innerTrack()->validFraction() < 0.8)
0903 return false;
0904 }
0905
0906 bool goodGlb = muon.isGlobalMuon() && muon.globalTrack()->normalizedChi2() < 3. &&
0907 muon.combinedQuality().chi2LocalPosition < 12. && muon.combinedQuality().trkKink < 20.;
0908
0909 return (segmentCompatibility(muon) > (goodGlb ? 0.303 : 0.451));
0910 }
0911
0912 bool muon::isSoftMuon(const reco::Muon& muon, const reco::Vertex& vtx, bool run2016_hip_mitigation) {
0913 bool muID = muon::isGoodMuon(muon, TMOneStationTight);
0914
0915 if (!muID)
0916 return false;
0917
0918 bool layers = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
0919 muon.innerTrack()->hitPattern().pixelLayersWithMeasurement() > 0;
0920
0921 bool ishighq = muon.innerTrack()->quality(reco::Track::highPurity);
0922
0923 bool ip =
0924 std::abs(muon.innerTrack()->dxy(vtx.position())) < 0.3 && std::abs(muon.innerTrack()->dz(vtx.position())) < 20.;
0925
0926 return layers && ip && (ishighq | run2016_hip_mitigation);
0927 }
0928
0929 bool muon::isHighPtMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
0930 if (!muon.isGlobalMuon())
0931 return false;
0932
0933 bool muValHits = (muon.globalTrack()->hitPattern().numberOfValidMuonHits() > 0 ||
0934 muon.tunePMuonBestTrack()->hitPattern().numberOfValidMuonHits() > 0);
0935
0936 bool muMatchedSt = muon.numberOfMatchedStations() > 1;
0937 if (!muMatchedSt) {
0938 if (muon.isTrackerMuon() && muon.numberOfMatchedStations() == 1) {
0939 if (muon.expectedNnumberOfMatchedStations() < 2 || !(muon.stationMask() == 1 || muon.stationMask() == 16) ||
0940 muon.numberOfMatchedRPCLayers() > 2)
0941 muMatchedSt = true;
0942 }
0943 }
0944
0945 bool muID = muValHits && muMatchedSt;
0946
0947 bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
0948 muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0;
0949
0950 bool momQuality = muon.tunePMuonBestTrack()->ptError() / muon.tunePMuonBestTrack()->pt() < 0.3;
0951
0952 bool ip =
0953 std::abs(muon.innerTrack()->dxy(vtx.position())) < 0.2 && std::abs(muon.innerTrack()->dz(vtx.position())) < 0.5;
0954
0955 return muID && hits && momQuality && ip;
0956 }
0957
0958 bool muon::isTrackerHighPtMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
0959 bool muID = muon.isTrackerMuon() && muon.track().isNonnull() && (muon.numberOfMatchedStations() > 1);
0960 if (!muID)
0961 return false;
0962
0963 bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
0964 muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0;
0965
0966 bool momQuality = muon.tunePMuonBestTrack()->ptError() < 0.3 * muon.tunePMuonBestTrack()->pt();
0967
0968 bool ip =
0969 std::abs(muon.innerTrack()->dxy(vtx.position())) < 0.2 && std::abs(muon.innerTrack()->dz(vtx.position())) < 0.5;
0970
0971 return muID && hits && momQuality && ip;
0972 }
0973
0974 int muon::sharedSegments(const reco::Muon& mu, const reco::Muon& mu2, unsigned int segmentArbitrationMask) {
0975 int ret = 0;
0976
0977
0978 for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = mu.matches().begin();
0979 chamberMatch != mu.matches().end();
0980 ++chamberMatch) {
0981 if (chamberMatch->segmentMatches.empty())
0982 continue;
0983 for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch2 = mu2.matches().begin();
0984 chamberMatch2 != mu2.matches().end();
0985 ++chamberMatch2) {
0986 if (chamberMatch2->segmentMatches.empty())
0987 continue;
0988 if (chamberMatch2->id() != chamberMatch->id())
0989 continue;
0990 for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();
0991 segmentMatch != chamberMatch->segmentMatches.end();
0992 ++segmentMatch) {
0993 if (!segmentMatch->isMask(segmentArbitrationMask))
0994 continue;
0995 for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch2 = chamberMatch2->segmentMatches.begin();
0996 segmentMatch2 != chamberMatch2->segmentMatches.end();
0997 ++segmentMatch2) {
0998 if (!segmentMatch2->isMask(segmentArbitrationMask))
0999 continue;
1000 if ((segmentMatch->cscSegmentRef.isNonnull() &&
1001 segmentMatch->cscSegmentRef == segmentMatch2->cscSegmentRef) ||
1002 (segmentMatch->dtSegmentRef.isNonnull() && segmentMatch->dtSegmentRef == segmentMatch2->dtSegmentRef)) {
1003 ++ret;
1004 }
1005 }
1006 }
1007 }
1008 }
1009
1010 return ret;
1011 }
1012
1013 bool outOfTimeMuon(const reco::Muon& muon) {
1014 const auto& combinedTime = muon.time();
1015 const auto& rpcTime = muon.rpcTime();
1016 bool combinedTimeIsOk = (combinedTime.nDof > 7);
1017 bool rpcTimeIsOk = (rpcTime.nDof > 1 && std::abs(rpcTime.timeAtIpInOutErr) < 0.001);
1018 bool outOfTime = false;
1019 if (rpcTimeIsOk) {
1020 if ((std::abs(rpcTime.timeAtIpInOut) > 10) && !(combinedTimeIsOk && std::abs(combinedTime.timeAtIpInOut) < 10))
1021 outOfTime = true;
1022 } else {
1023 if (combinedTimeIsOk && (combinedTime.timeAtIpInOut > 20 || combinedTime.timeAtIpInOut < -45))
1024 outOfTime = true;
1025 }
1026 return outOfTime;
1027 }
1028
1029 reco::Muon::Selector muon::makeSelectorBitset(reco::Muon const& muon,
1030 reco::Vertex const* vertex,
1031 bool run2016_hip_mitigation) {
1032
1033 unsigned int selectors = muon.selectors();
1034
1035 double chIso = muon.pfIsolationR04().sumChargedHadronPt;
1036 double nIso = muon.pfIsolationR04().sumNeutralHadronEt;
1037 double phoIso = muon.pfIsolationR04().sumPhotonEt;
1038 double puIso = muon.pfIsolationR04().sumPUPt;
1039 double dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5 * puIso, 0.);
1040 double dbCorrectedRelIso = dbCorrectedIsolation / muon.pt();
1041 double tkRelIso = muon.isolationR03().sumPt / muon.pt();
1042
1043
1044 if (muon::isLooseMuon(muon))
1045 selectors |= (1UL << reco::Muon::CutBasedIdLoose);
1046 if (vertex) {
1047 if (muon::isTightMuon(muon, *vertex))
1048 selectors |= (1UL << reco::Muon::CutBasedIdTight);
1049 if (muon::isSoftMuon(muon, *vertex, run2016_hip_mitigation))
1050 selectors |= (1UL << reco::Muon::SoftCutBasedId);
1051 if (muon::isHighPtMuon(muon, *vertex))
1052 selectors |= (1UL << reco::Muon::CutBasedIdGlobalHighPt);
1053 if (muon::isTrackerHighPtMuon(muon, *vertex))
1054 selectors |= (1UL << reco::Muon::CutBasedIdTrkHighPt);
1055 }
1056 if (muon::isMediumMuon(muon, run2016_hip_mitigation)) {
1057 selectors |= (1UL << reco::Muon::CutBasedIdMedium);
1058 if (vertex and std::abs(muon.muonBestTrack()->dz(vertex->position())) < 0.1 and
1059 std::abs(muon.muonBestTrack()->dxy(vertex->position())) < 0.02)
1060 selectors |= (1UL << reco::Muon::CutBasedIdMediumPrompt);
1061 }
1062
1063
1064 if (dbCorrectedRelIso < 0.40)
1065 selectors |= (1UL << reco::Muon::PFIsoVeryLoose);
1066 if (dbCorrectedRelIso < 0.25)
1067 selectors |= (1UL << reco::Muon::PFIsoLoose);
1068 if (dbCorrectedRelIso < 0.20)
1069 selectors |= (1UL << reco::Muon::PFIsoMedium);
1070 if (dbCorrectedRelIso < 0.15)
1071 selectors |= (1UL << reco::Muon::PFIsoTight);
1072 if (dbCorrectedRelIso < 0.10)
1073 selectors |= (1UL << reco::Muon::PFIsoVeryTight);
1074 if (dbCorrectedRelIso < 0.05)
1075 selectors |= (1UL << reco::Muon::PFIsoVeryVeryTight);
1076
1077
1078 if (tkRelIso < 0.10)
1079 selectors |= (1UL << reco::Muon::TkIsoLoose);
1080 if (tkRelIso < 0.05)
1081 selectors |= (1UL << reco::Muon::TkIsoTight);
1082
1083
1084 if (isLooseTriggerMuon(muon))
1085 selectors |= (1UL << reco::Muon::TriggerIdLoose);
1086
1087
1088 if (!outOfTimeMuon(muon))
1089 selectors |= (1UL << reco::Muon::InTimeMuon);
1090
1091 return static_cast<reco::Muon::Selector>(selectors);
1092 }