MuonReco/src/MuonSelectors.cc

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     // in case of unrecognized selection type
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     // in case of unrecognized selection type
0034     if (!found)
0035       throw cms::Exception("MuonSelectorError") << label << " is not a recognized reco::Muon::Selector";
0036     return value;
0037   }
0038 }  // namespace muon
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 // ------------ method to calculate the calo compatibility for a track with matched muon info  ------------
0058 float muon::caloCompatibility(const reco::Muon& muon) { return muon.caloCompatibility(); }
0059
0060 // ------------ method to calculate the segment compatibility for a track with matched muon info  ------------
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     // *** fill local info for this muon (do some counting) ***;
0078     // ************** begin ***********************************;
0079     if (i <= 4) {  // this is the section for the DTs
0080       float thisTrackDist = muon.trackDist(i, 1, arbitrationType);
0081       if (thisTrackDist < 999999) {  //current "raw" info that a track is close to a chamber
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       //current "raw" info that a segment is matched to the current track
0090       if (muon.segmentX(i, 1, arbitrationType) < 999999) {
0091         ++nr_of_stations_with_segment;
0092         station_has_segmentmatch[i - 1] = 1;
0093       }
0094     } else {  // this is the section for the CSCs
0095       float thisTrackDist = muon.trackDist(i - 4, 2, arbitrationType);
0096       if (thisTrackDist < 999999) {  //current "raw" info that a track is close to a chamber
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       //current "raw" info that a segment is matched to the current track
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     // rough estimation of chamber border efficiency (should be parametrized better, this is just a quick guess):
0111     // TF1 * merf = new TF1("merf","-0.5*(TMath::Erf(x/6.)-1)",-100,100);
0112     // use above value to "unpunish" missing segment if close to border, i.e. rather than not adding any weight, add
0113     // the one from the function. Only for dist ~> -10 cm, else full punish!.
0114
0115     // ********************************************************;
0116     // *** fill local info for this muon (do some counting) ***;
0117     // ************** end *************************************;
0118   }
0119
0120   // ********************************************************;
0121   // *** calculate weights for each station *****************;
0122   // ************** begin ***********************************;
0123   //    const float slope = 0.5;
0124   //    const float attenuate_weight_regain = 1.;
0125   // if attenuate_weight_regain < 1., additional punishment if track is close to boundary and no segment
0126   const float attenuate_weight_regain = 0.5;
0127
0128   for (int i = 1; i <= 8; ++i) {  // loop over all possible stations
0129
0130     // first set all weights if a station has been crossed
0131     // later penalize if a station did not have a matching segment
0132
0133     //old logic      if(station_has_segmentmatch[i-1] > 0 ) { // the track has an associated segment at the current station
0134     if (station_was_crossed[i - 1] > 0) {  // the track crossed this chamber (or was nearby)
0135       // - Apply a weight depending on the "depth" of the muon passage.
0136       // - The station_weight is later reduced for stations with badly matched segments.
0137       // - Even if there is no segment but the track passes close to a chamber boundary, the
0138       //   weight is set non zero and can go up to 0.5 of the full weight if the track is quite
0139       //   far from any station.
0140       ++position_in_stations;
0141
0142       switch (nr_of_stations_crossed) {  // define different weights depending on how many stations were 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           //    LogTrace("MuonIdentification")<<"            // Message: A muon candidate track has more than 4 stations with matching segments.";
0173           //    LogTrace("MuonIdentification")<<"            // Did not expect this - please let me know: ibloch@fnal.gov";
0174           // for all other cases
0175           station_weight[i - 1] = 1.f / nr_of_stations_crossed;
0176       }
0177
0178       if (use_weight_regain_at_chamber_boundary) {  // reconstitute some weight if there is no match but the segment is close to a boundary:
0179         if (station_has_segmentmatch[i - 1] <= 0 && stations_w_track_at_boundary[i - 1] != 0.) {
0180           // if segment is not present but track in inefficient region, do not count as "missing match" but add some reduced weight.
0181           // original "match weight" is currently reduced by at least attenuate_weight_regain, variing with an error function down to 0 if the track is
0182           // inside the chamber.
0183           // remark: the additional scale of 0.5 normalizes Err to run from 0 to 1 in y
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) {  // no segment match and track well inside chamber
0188           // full penalization
0189           station_weight[i - 1] = 0.f;
0190         }
0191       } else {  // always fully penalize tracks with no matching segment, whether the segment is close to the boundary or not.
0192         if (station_has_segmentmatch[i - 1] <= 0)
0193           station_weight[i - 1] = 0.f;
0194       }
0195
0196       // if track has matching segment, but the matching is not high quality, penalize
0197       if (station_has_segmentmatch[i - 1] > 0) {
0198         if (i <= 4) {  // we are in the DTs
0199           if (muon.dY(i, 1, arbitrationType) < 999999.f &&
0200               muon.dX(i, 1, arbitrationType) < 999999.f) {  // have both X and Y match
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               // reduce weight
0207               if (use_match_dist_penalty) {
0208                 // only use pull if 3 sigma is not smaller than 3 cm
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) {  // has no match in Y
0218             // has a match in X. Pull larger that 1 to avoid increasing the weight (just penalize, don't anti-penalize)
0219             if (muon.pullX(i, 1, arbitrationType) > 1.f) {
0220               // reduce weight
0221               if (use_match_dist_penalty) {
0222                 // only use pull if 3 sigma is not smaller than 3 cm
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 {  // has no match in X
0232             // has a match in Y. Pull larger that 1 to avoid increasing the weight (just penalize, don't anti-penalize)
0233             if (muon.pullY(i, 1, arbitrationType) > 1.f) {
0234               // reduce weight
0235               if (use_match_dist_penalty) {
0236                 // only use pull if 3 sigma is not smaller than 3 cm
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 {  // We are in the CSCs
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             // reduce weight
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               // only use pull if 3 sigma is not smaller than 3 cm
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       // Thoughts:
0267       // - should penalize if the segment has only x OR y info
0268       // - should also use the segment direction, as it now works!
0269
0270     } else {  // track did not pass a chamber in this station - just reset weight
0271       station_weight[i - 1] = 0.f;
0272     }
0273
0274     //increment final weight for muon:
0275     full_weight += station_weight[i - 1];
0276   }
0277
0278   // if we don't expect any matches, we set the compatibility to
0279   // 0.5 as the track is as compatible with a muon as it is with
0280   // background - we should maybe rather set it to -0.5!
0281   if (nr_of_stations_crossed == 0) {
0282     //      full_weight = attenuate_weight_regain*0.5;
0283     full_weight = 0.5f;
0284   }
0285
0286   // ********************************************************;
0287   // *** calculate weights for each station *****************;
0288   // ************** end *************************************;
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       // Simplistic first cut in the 2D segment- vs calo-compatibility plane. Will have to be refined!
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       //    LogTrace("MuonIdentification")<<"            // Invalid Algorithm Type called!";
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     // To satisfy my own paranoia, if the user specifies that the
0336     // minimum number of matches is zero, then return true.
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     // Require that there be at least a minimum number of segments
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     // Make sure the minimum number of matches is not greater than
0355     // the number of required stations but still greater than zero
0356     if (syncMinNMatchesNRequiredStationsInBarrelOnly) {
0357       // Note that we only do this in the barrel region!
0358       if (std::abs(muon.eta()) < 1.2) {
0359         if (minNumberOfMatches > numRequiredStations)
0360           minNumberOfMatches = numRequiredStations;
0361         if (minNumberOfMatches < 1)  //SK: this only happens for negative values
0362           minNumberOfMatches = 1;
0363       }
0364     } else {
0365       if (minNumberOfMatches > numRequiredStations)
0366         minNumberOfMatches = numRequiredStations;
0367       if (minNumberOfMatches < 1)  //SK: this only happens for negative values
0368         minNumberOfMatches = 1;
0369     }
0370
0371     if (numSegs >= minNumberOfMatches)
0372       goodMuon = true;
0373
0374     // Require that last required station have segment
0375     // If there are zero required stations keep track
0376     // of the last station with a segment so that we may
0377     // apply the quality cuts below to it instead
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     // Impose pull cuts on last segment
0403     int station = 0, detector = 0;
0404     station = lastSegBit < 4 ? lastSegBit + 1 : lastSegBit - 3;
0405     detector = lastSegBit < 4 ? 1 : 2;
0406
0407     // Check x information
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     // Is this a tight algorithm, i.e. do we bother to check y information?
0416     if (maxAbsDy < 999999) {  // really if maxAbsDy < 1E9 as currently defined
0417
0418       // Check y information
0419       if (detector == 2) {  // CSC
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         // In DT, if this is a "Tight" algorithm and the last segment is
0429         // missing y information (always the case in station 4!!!), impose
0430         // respective cuts on the next station in the stationMask that has
0431         // a segment with y information.  If there are no segments with y
0432         // information then there is nothing to penalize. Should we
0433         // penalize in Tight for having zero segments with y information?
0434         // That is the fundamental question.  Of course I am being uber
0435         // paranoid; if this is a good muon then there will probably be at
0436         // least one segment with y information but not always.  Suppose
0437         // somehow a muon only creates segments in station 4, then we
0438         // definitely do not want to require that there be at least one
0439         // segment with y information because we will lose it completely.
0440         //
0441
0442         for (int stationIdx = station; stationIdx > 0; --stationIdx) {
0443           if (!(theStationMask & 1 << (stationIdx - 1)))  // don't bother if the station is not in the stationMask
0444             continue;
0445
0446           if (muon.dY(stationIdx, 1, arbitrationType) > 999998)  // no y-information
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           // If we get this far then great this is a good muon
0458           return true;
0459         }
0460       }
0461     }
0462
0463     return goodMuon;
0464   }  // TMLastStation
0465
0466   // TMOneStation requires only that there be one "good" segment, regardless
0467   // of the required stations.  We do not penalize if there are absolutely zero
0468   // segments with y information in the Tight algorithm.  Maybe I'm being
0469   // paranoid but so be it.  If it's really a good muon then we will probably
0470   // find at least one segment with both x and y information but you never
0471   // know, and I don't want to deal with a potential inefficiency in the DT
0472   // like we did with the original TMLastStation.  Incidentally, not penalizing
0473   // for total lack of y information in the Tight algorithm is what is done in
0474   // the new TMLastStation
0475   //
0476   if (type == TMOneStation) {
0477     unsigned int theStationMask = muon.stationMask(arbitrationType);
0478
0479     // Of course there must be at least one segment
0480     if (!theStationMask)
0481       return false;
0482
0483     int station = 0, detector = 0;
0484     // Keep track of whether or not there is a DT segment with y information.
0485     // In the end, if it turns out there are absolutely zero DT segments with
0486     // y information, require only that there was a segment with good x info.
0487     // This of course only applies to the Tight algorithms.
0488     bool existsGoodDTSegX = false;
0489     bool existsDTSegY = false;
0490
0491     // Impose cuts on the segments in the station mask until we find a good one
0492     // Might as well start with the lowest bit to speed things up.
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         // Is this a tight algorithm?  If yes, use y information
0506         if (maxAbsDy < 999999) {
0507           if (detector == 2) {  // CSC
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)  // no y-information
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         // If we get this far then great this is a good muon
0527         return true;
0528       }
0529
0530     // If we get this far then for sure there are no "good" CSC segments. For
0531     // DT, check if there were any segments with y information.  If there
0532     // were none, but there was a segment with good x, then we're happy. If
0533     // there WERE segments with y information, then they must have been shit
0534     // since we are here so fail it.  Of course, if this is a Loose algorithm
0535     // then fail immediately since if we had good x we would already have
0536     // returned true
0537     if (maxAbsDy < 999999) {
0538       if (existsDTSegY)
0539         return false;
0540       else if (existsGoodDTSegX)
0541         return true;
0542     } else
0543       return false;
0544   }  // TMOneStation
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   }  // RPCMu
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   }  // ME0Mu
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   }  // GEMMu
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       // For "Loose" algorithms we choose maximum y quantity cuts of 1E9 instead of
0674       // 9999 as before.  We do this because the muon methods return 999999 (note
0675       // there are six 9's) when the requested information is not available.  For
0676       // example, if a muon fails to traverse the z measuring superlayer in a station
0677       // in the DT, then all methods involving segmentY in this station return
0678       // 999999 to demonstrate that the information is missing.  In order to not
0679       // penalize muons for missing y information in Loose algorithms where we do
0680       // not care at all about y information, we raise these limits.  In the
0681       // TMLastStation and TMOneStation algorithms we actually use this huge number
0682       // to determine whether to consider y information at all.
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       //compatibility loose
0716     case muon::TM2DCompatibilityLoose:
0717       return muon.isTrackerMuon() && isGoodMuon(muon, TM2DCompatibility, 0.7, arbitrationType);
0718       break;
0719       //compatibility tight
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       // if ( (chamber1->segmentMatches.empty() || chamber2->segmentMatches.empty()) ) continue;
0803
0804       // handle case where both muons have information about the same chamber
0805       // here we know how close they are
0806       if (chamber1->id == chamber2->id) {
0807         // found the same chamber
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         // check if tracks are pointing into overlaping region of the CSC detector
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         // FIXME: we don't handle 18->1; 36->1 transitions since
0844         // I don't know how to check for sure how many chambers
0845         // are there. Probably need to hard code some checks.
0846
0847         // Now we have to make sure that both tracks are close to an edge
0848         // FIXME: ignored Y coordinate for now
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) {  // check if the same edge
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   // Requirements:
0868   // - no depencence on information not availabe in the muon object
0869   // - use only robust inputs
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   // Will do with a stupid double loop, since creating and filling a map is probably _more_ inefficient for a single lookup.
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           }  // is the same
1009         }    // segment of mu2 in chamber
1010       }      // segment of mu1 in chamber
1011     }        // chamber of mu2
1012   }          // chamber of mu1
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   // https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideMuonIdRun2
1037   unsigned int selectors = muon.selectors();
1038   // Compute Id and Isolation variables
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   // Base selectors
1048   if (muon::isLooseMuon(muon))
1049     selectors |= reco::Muon::CutBasedIdLoose;
1050   if (vertex) {
1051     if (muon::isTightMuon(muon, *vertex))
1052       selectors |= reco::Muon::CutBasedIdTight;
1053     if (muon::isSoftMuon(muon, *vertex, run2016_hip_mitigation))
1054       selectors |= reco::Muon::SoftCutBasedId;
1055     if (muon::isHighPtMuon(muon, *vertex))
1056       selectors |= reco::Muon::CutBasedIdGlobalHighPt;
1057     if (muon::isTrackerHighPtMuon(muon, *vertex))
1058       selectors |= reco::Muon::CutBasedIdTrkHighPt;
1059   }
1060   if (muon::isMediumMuon(muon, run2016_hip_mitigation)) {
1061     selectors |= 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 |= reco::Muon::CutBasedIdMediumPrompt;
1065   }
1066
1067   // PF isolation
1068   if (dbCorrectedRelIso < 0.40)
1069     selectors |= reco::Muon::PFIsoVeryLoose;
1070   if (dbCorrectedRelIso < 0.25)
1071     selectors |= reco::Muon::PFIsoLoose;
1072   if (dbCorrectedRelIso < 0.20)
1073     selectors |= reco::Muon::PFIsoMedium;
1074   if (dbCorrectedRelIso < 0.15)
1075     selectors |= reco::Muon::PFIsoTight;
1076   if (dbCorrectedRelIso < 0.10)
1077     selectors |= reco::Muon::PFIsoVeryTight;
1078   if (dbCorrectedRelIso < 0.05)
1079     selectors |= reco::Muon::PFIsoVeryVeryTight;
1080
1081   // Tracker isolation
1082   if (tkRelIso < 0.10)
1083     selectors |= reco::Muon::TkIsoLoose;
1084   if (tkRelIso < 0.05)
1085     selectors |= reco::Muon::TkIsoTight;
1086
1087   // Trigger selectors
1088   if (isLooseTriggerMuon(muon))
1089     selectors |= reco::Muon::TriggerIdLoose;
1090
1091   // Timing
1092   if (!outOfTimeMuon(muon))
1093     selectors |= reco::Muon::InTimeMuon;
1094
1095   return static_cast<reco::Muon::Selector>(selectors);
1096 }