Back to home page

Project CMSSW displayed by LXR

 
 

    


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     // 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   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     // *** fill local info for this muon (do some counting) ***;
0077     // ************** begin ***********************************;
0078     if (i <= 4) {  // this is the section for the DTs
0079       float thisTrackDist = muon.trackDist(i, 1, arbitrationType);
0080       if (thisTrackDist < 999999) {  //current "raw" info that a track is close to a chamber
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       //current "raw" info that a segment is matched to the current track
0089       if (muon.segmentX(i, 1, arbitrationType) < 999999) {
0090         station_has_segmentmatch[i - 1] = 1;
0091       }
0092     } else {  // this is the section for the CSCs
0093       float thisTrackDist = muon.trackDist(i - 4, 2, arbitrationType);
0094       if (thisTrackDist < 999999) {  //current "raw" info that a track is close to a chamber
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       //current "raw" info that a segment is matched to the current track
0103       if (muon.segmentX(i - 4, 2, arbitrationType) < 999999) {
0104         station_has_segmentmatch[i - 1] = 1;
0105       }
0106     }
0107     // rough estimation of chamber border efficiency (should be parametrized better, this is just a quick guess):
0108     // TF1 * merf = new TF1("merf","-0.5*(TMath::Erf(x/6.)-1)",-100,100);
0109     // use above value to "unpunish" missing segment if close to border, i.e. rather than not adding any weight, add
0110     // the one from the function. Only for dist ~> -10 cm, else full punish!.
0111 
0112     // ********************************************************;
0113     // *** fill local info for this muon (do some counting) ***;
0114     // ************** end *************************************;
0115   }
0116 
0117   // ********************************************************;
0118   // *** calculate weights for each station *****************;
0119   // ************** begin ***********************************;
0120   //    const float slope = 0.5;
0121   //    const float attenuate_weight_regain = 1.;
0122   // if attenuate_weight_regain < 1., additional punishment if track is close to boundary and no segment
0123   const float attenuate_weight_regain = 0.5;
0124 
0125   for (int i = 1; i <= 8; ++i) {  // loop over all possible stations
0126 
0127     // first set all weights if a station has been crossed
0128     // later penalize if a station did not have a matching segment
0129 
0130     //old logic      if(station_has_segmentmatch[i-1] > 0 ) { // the track has an associated segment at the current station
0131     if (station_was_crossed[i - 1] > 0) {  // the track crossed this chamber (or was nearby)
0132       // - Apply a weight depending on the "depth" of the muon passage.
0133       // - The station_weight is later reduced for stations with badly matched segments.
0134       // - Even if there is no segment but the track passes close to a chamber boundary, the
0135       //   weight is set non zero and can go up to 0.5 of the full weight if the track is quite
0136       //   far from any station.
0137       ++position_in_stations;
0138 
0139       switch (nr_of_stations_crossed) {  // define different weights depending on how many stations were 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           //    LogTrace("MuonIdentification")<<"            // Message: A muon candidate track has more than 4 stations with matching segments.";
0170           //    LogTrace("MuonIdentification")<<"            // Did not expect this - please let me know: ibloch@fnal.gov";
0171           // for all other cases
0172           station_weight[i - 1] = 1.f / nr_of_stations_crossed;
0173       }
0174 
0175       if (use_weight_regain_at_chamber_boundary) {  // reconstitute some weight if there is no match but the segment is close to a boundary:
0176         if (station_has_segmentmatch[i - 1] <= 0 && stations_w_track_at_boundary[i - 1] != 0.) {
0177           // if segment is not present but track in inefficient region, do not count as "missing match" but add some reduced weight.
0178           // original "match weight" is currently reduced by at least attenuate_weight_regain, variing with an error function down to 0 if the track is
0179           // inside the chamber.
0180           // remark: the additional scale of 0.5 normalizes Err to run from 0 to 1 in y
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) {  // no segment match and track well inside chamber
0185           // full penalization
0186           station_weight[i - 1] = 0.f;
0187         }
0188       } else {  // always fully penalize tracks with no matching segment, whether the segment is close to the boundary or not.
0189         if (station_has_segmentmatch[i - 1] <= 0)
0190           station_weight[i - 1] = 0.f;
0191       }
0192 
0193       // if track has matching segment, but the matching is not high quality, penalize
0194       if (station_has_segmentmatch[i - 1] > 0) {
0195         if (i <= 4) {  // we are in the DTs
0196           if (muon.dY(i, 1, arbitrationType) < 999999.f &&
0197               muon.dX(i, 1, arbitrationType) < 999999.f) {  // have both X and Y match
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               // reduce weight
0204               if (use_match_dist_penalty) {
0205                 // only use pull if 3 sigma is not smaller than 3 cm
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) {  // has no match in Y
0215             // has a match in X. Pull larger that 1 to avoid increasing the weight (just penalize, don't anti-penalize)
0216             if (muon.pullX(i, 1, arbitrationType) > 1.f) {
0217               // reduce weight
0218               if (use_match_dist_penalty) {
0219                 // only use pull if 3 sigma is not smaller than 3 cm
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 {  // has no match in X
0229             // has a match in Y. Pull larger that 1 to avoid increasing the weight (just penalize, don't anti-penalize)
0230             if (muon.pullY(i, 1, arbitrationType) > 1.f) {
0231               // reduce weight
0232               if (use_match_dist_penalty) {
0233                 // only use pull if 3 sigma is not smaller than 3 cm
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 {  // We are in the CSCs
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             // reduce weight
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               // only use pull if 3 sigma is not smaller than 3 cm
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       // Thoughts:
0264       // - should penalize if the segment has only x OR y info
0265       // - should also use the segment direction, as it now works!
0266 
0267     } else {  // track did not pass a chamber in this station - just reset weight
0268       station_weight[i - 1] = 0.f;
0269     }
0270 
0271     //increment final weight for muon:
0272     full_weight += station_weight[i - 1];
0273   }
0274 
0275   // if we don't expect any matches, we set the compatibility to
0276   // 0.5 as the track is as compatible with a muon as it is with
0277   // background - we should maybe rather set it to -0.5!
0278   if (nr_of_stations_crossed == 0) {
0279     //      full_weight = attenuate_weight_regain*0.5;
0280     full_weight = 0.5f;
0281   }
0282 
0283   // ********************************************************;
0284   // *** calculate weights for each station *****************;
0285   // ************** end *************************************;
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       // Simplistic first cut in the 2D segment- vs calo-compatibility plane. Will have to be refined!
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       //    LogTrace("MuonIdentification")<<"            // Invalid Algorithm Type called!";
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     // To satisfy my own paranoia, if the user specifies that the
0333     // minimum number of matches is zero, then return true.
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     // Require that there be at least a minimum number of segments
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     // Make sure the minimum number of matches is not greater than
0352     // the number of required stations but still greater than zero
0353     if (syncMinNMatchesNRequiredStationsInBarrelOnly) {
0354       // Note that we only do this in the barrel region!
0355       if (std::abs(muon.eta()) < 1.2) {
0356         if (minNumberOfMatches > numRequiredStations)
0357           minNumberOfMatches = numRequiredStations;
0358         if (minNumberOfMatches < 1)  //SK: this only happens for negative values
0359           minNumberOfMatches = 1;
0360       }
0361     } else {
0362       if (minNumberOfMatches > numRequiredStations)
0363         minNumberOfMatches = numRequiredStations;
0364       if (minNumberOfMatches < 1)  //SK: this only happens for negative values
0365         minNumberOfMatches = 1;
0366     }
0367 
0368     if (numSegs >= minNumberOfMatches)
0369       goodMuon = true;
0370 
0371     // Require that last required station have segment
0372     // If there are zero required stations keep track
0373     // of the last station with a segment so that we may
0374     // apply the quality cuts below to it instead
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     // Impose pull cuts on last segment
0400     int station = 0, detector = 0;
0401     station = lastSegBit < 4 ? lastSegBit + 1 : lastSegBit - 3;
0402     detector = lastSegBit < 4 ? 1 : 2;
0403 
0404     // Check x information
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     // Is this a tight algorithm, i.e. do we bother to check y information?
0413     if (maxAbsDy < 999999) {  // really if maxAbsDy < 1E9 as currently defined
0414 
0415       // Check y information
0416       if (detector == 2) {  // CSC
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         // In DT, if this is a "Tight" algorithm and the last segment is
0426         // missing y information (always the case in station 4!!!), impose
0427         // respective cuts on the next station in the stationMask that has
0428         // a segment with y information.  If there are no segments with y
0429         // information then there is nothing to penalize. Should we
0430         // penalize in Tight for having zero segments with y information?
0431         // That is the fundamental question.  Of course I am being uber
0432         // paranoid; if this is a good muon then there will probably be at
0433         // least one segment with y information but not always.  Suppose
0434         // somehow a muon only creates segments in station 4, then we
0435         // definitely do not want to require that there be at least one
0436         // segment with y information because we will lose it completely.
0437         //
0438 
0439         for (int stationIdx = station; stationIdx > 0; --stationIdx) {
0440           if (!(theStationMask & 1 << (stationIdx - 1)))  // don't bother if the station is not in the stationMask
0441             continue;
0442 
0443           if (muon.dY(stationIdx, 1, arbitrationType) > 999998)  // no y-information
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           // If we get this far then great this is a good muon
0455           return true;
0456         }
0457       }
0458     }
0459 
0460     return goodMuon;
0461   }  // TMLastStation
0462 
0463   // TMOneStation requires only that there be one "good" segment, regardless
0464   // of the required stations.  We do not penalize if there are absolutely zero
0465   // segments with y information in the Tight algorithm.  Maybe I'm being
0466   // paranoid but so be it.  If it's really a good muon then we will probably
0467   // find at least one segment with both x and y information but you never
0468   // know, and I don't want to deal with a potential inefficiency in the DT
0469   // like we did with the original TMLastStation.  Incidentally, not penalizing
0470   // for total lack of y information in the Tight algorithm is what is done in
0471   // the new TMLastStation
0472   //
0473   if (type == TMOneStation) {
0474     unsigned int theStationMask = muon.stationMask(arbitrationType);
0475 
0476     // Of course there must be at least one segment
0477     if (!theStationMask)
0478       return false;
0479 
0480     int station = 0, detector = 0;
0481     // Keep track of whether or not there is a DT segment with y information.
0482     // In the end, if it turns out there are absolutely zero DT segments with
0483     // y information, require only that there was a segment with good x info.
0484     // This of course only applies to the Tight algorithms.
0485     bool existsGoodDTSegX = false;
0486     bool existsDTSegY = false;
0487 
0488     // Impose cuts on the segments in the station mask until we find a good one
0489     // Might as well start with the lowest bit to speed things up.
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         // Is this a tight algorithm?  If yes, use y information
0503         if (maxAbsDy < 999999) {
0504           if (detector == 2) {  // CSC
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)  // no y-information
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         // If we get this far then great this is a good muon
0524         return true;
0525       }
0526 
0527     // If we get this far then for sure there are no "good" CSC segments. For
0528     // DT, check if there were any segments with y information.  If there
0529     // were none, but there was a segment with good x, then we're happy. If
0530     // there WERE segments with y information, then they must have been shit
0531     // since we are here so fail it.  Of course, if this is a Loose algorithm
0532     // then fail immediately since if we had good x we would already have
0533     // returned true
0534     if (maxAbsDy < 999999) {
0535       if (existsDTSegY)
0536         return false;
0537       else if (existsGoodDTSegX)
0538         return true;
0539     } else
0540       return false;
0541   }  // TMOneStation
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   }  // RPCMu
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   }  // ME0Mu
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   }  // GEMMu
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       // For "Loose" algorithms we choose maximum y quantity cuts of 1E9 instead of
0671       // 9999 as before.  We do this because the muon methods return 999999 (note
0672       // there are six 9's) when the requested information is not available.  For
0673       // example, if a muon fails to traverse the z measuring superlayer in a station
0674       // in the DT, then all methods involving segmentY in this station return
0675       // 999999 to demonstrate that the information is missing.  In order to not
0676       // penalize muons for missing y information in Loose algorithms where we do
0677       // not care at all about y information, we raise these limits.  In the
0678       // TMLastStation and TMOneStation algorithms we actually use this huge number
0679       // to determine whether to consider y information at all.
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       //compatibility loose
0713     case muon::TM2DCompatibilityLoose:
0714       return muon.isTrackerMuon() && isGoodMuon(muon, TM2DCompatibility, 0.7, arbitrationType);
0715       break;
0716       //compatibility tight
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       // if ( (chamber1->segmentMatches.empty() || chamber2->segmentMatches.empty()) ) continue;
0800 
0801       // handle case where both muons have information about the same chamber
0802       // here we know how close they are
0803       if (chamber1->id == chamber2->id) {
0804         // found the same chamber
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         // check if tracks are pointing into overlaping region of the CSC detector
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         // FIXME: we don't handle 18->1; 36->1 transitions since
0841         // I don't know how to check for sure how many chambers
0842         // are there. Probably need to hard code some checks.
0843 
0844         // Now we have to make sure that both tracks are close to an edge
0845         // FIXME: ignored Y coordinate for now
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) {  // check if the same edge
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   // Requirements:
0865   // - no depencence on information not availabe in the muon object
0866   // - use only robust inputs
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   // Will do with a stupid double loop, since creating and filling a map is probably _more_ inefficient for a single lookup.
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           }  // is the same
1005         }  // segment of mu2 in chamber
1006       }  // segment of mu1 in chamber
1007     }  // chamber of mu2
1008   }  // chamber of mu1
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   // https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideMuonIdRun2
1033   unsigned int selectors = muon.selectors();
1034   // Compute Id and Isolation variables
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   // Base selectors
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   // PF isolation
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   // Tracker isolation
1078   if (tkRelIso < 0.10)
1079     selectors |= (1UL << reco::Muon::TkIsoLoose);
1080   if (tkRelIso < 0.05)
1081     selectors |= (1UL << reco::Muon::TkIsoTight);
1082 
1083   // Trigger selectors
1084   if (isLooseTriggerMuon(muon))
1085     selectors |= (1UL << reco::Muon::TriggerIdLoose);
1086 
1087   // Timing
1088   if (!outOfTimeMuon(muon))
1089     selectors |= (1UL << reco::Muon::InTimeMuon);
1090 
1091   return static_cast<reco::Muon::Selector>(selectors);
1092 }