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