Project CMSSW displayed by LXR CMSSW_14_0_X_2023-12-10-2300/​DataFormats/​MuonReco/​src/​MuonSelectors.cc	Source navigation Diff markup Identifier search General search
	Version: CMSSW_14_0_X_2023-12-10-2300 CMSSW_14_0_X_2023-12-08-2300 CMSSW_14_0_X_2023-12-07-2300 CMSSW_14_0_X_2023-12-06-2300 CMSSW_14_0_X_2023-12-05-2300 CMSSW_14_0_X_2023-12-04-2300 Revert   Change

File indexing completed on 2023-10-25 09:39:35

 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
 #include "DataFormats/TrackReco/interface/Track.h"
 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
 #include "DataFormats/VertexReco/interface/Vertex.h"
 #include "DataFormats/MuonReco/interface/MuonRPCHitMatch.h"
 
 namespace muon {
   SelectionType selectionTypeFromString(const std::string& label) {
     SelectionType value = (SelectionType)-1;
     bool found = false;
     for (int i = 0; selectionTypeStringToEnumMap[i].label && (!found); ++i)
       if (!strcmp(label.c_str(), selectionTypeStringToEnumMap[i].label)) {
         found = true;
         value = selectionTypeStringToEnumMap[i].value;
       }
 
     // in case of unrecognized selection type
     if (!found)
       throw cms::Exception("MuonSelectorError") << label << " is not a recognized SelectionType";
     return value;
   }
 
   reco::Muon::Selector selectorFromString(const std::string& label) {
     reco::Muon::Selector value = (reco::Muon::Selector)-1;
     bool found = false;
     for (int i = 0; selectorStringToEnumMap[i].label && (!found); ++i)
       if (!strcmp(label.c_str(), selectorStringToEnumMap[i].label)) {
         found = true;
         value = selectorStringToEnumMap[i].value;
       }
 
     // in case of unrecognized selection type
     if (!found)
       throw cms::Exception("MuonSelectorError") << label << " is not a recognized reco::Muon::Selector";
     return value;
   }
 }  // namespace muon
 
 unsigned int muon::RequiredStationMask(const reco::Muon& muon,
                                        double maxChamberDist,
                                        double maxChamberDistPull,
                                        reco::Muon::ArbitrationType arbitrationType) {
   unsigned int theMask = 0;
 
   for (int stationIdx = 1; stationIdx < 5; ++stationIdx) {
     for (int detectorIdx = 1; detectorIdx < 3; ++detectorIdx) {
       float dist = muon.trackDist(stationIdx, detectorIdx, arbitrationType);
       if (dist < maxChamberDist &&
           dist < maxChamberDistPull * muon.trackDistErr(stationIdx, detectorIdx, arbitrationType))
         theMask += 1 << ((stationIdx - 1) + 4 * (detectorIdx - 1));
     }
   }
   return theMask;
 }
 
 // ------------ method to calculate the calo compatibility for a track with matched muon info  ------------
 float muon::caloCompatibility(const reco::Muon& muon) { return muon.caloCompatibility(); }
 
 // ------------ method to calculate the segment compatibility for a track with matched muon info  ------------
 float muon::segmentCompatibility(const reco::Muon& muon, reco::Muon::ArbitrationType arbitrationType) {
   bool use_weight_regain_at_chamber_boundary = true;
   bool use_match_dist_penalty = true;
 
   int nr_of_stations_crossed = 0;
   std::vector<int> stations_w_track(8);
   std::vector<int> station_has_segmentmatch(8);
   std::vector<int> station_was_crossed(8);
   std::vector<float> stations_w_track_at_boundary(8);
   std::vector<float> station_weight(8);
   int position_in_stations = 0;
   float full_weight = 0.;
 
   for (int i = 1; i <= 8; ++i) {
     // ********************************************************;
     // *** fill local info for this muon (do some counting) ***;
     // ************** begin ***********************************;
     if (i <= 4) {  // this is the section for the DTs
       float thisTrackDist = muon.trackDist(i, 1, arbitrationType);
       if (thisTrackDist < 999999) {  //current "raw" info that a track is close to a chamber
         ++nr_of_stations_crossed;
         station_was_crossed[i - 1] = 1;
         if (thisTrackDist > -10.)
           stations_w_track_at_boundary[i - 1] = thisTrackDist;
         else
           stations_w_track_at_boundary[i - 1] = 0.;
       }
       //current "raw" info that a segment is matched to the current track
       if (muon.segmentX(i, 1, arbitrationType) < 999999) {
         station_has_segmentmatch[i - 1] = 1;
       }
     } else {  // this is the section for the CSCs
       float thisTrackDist = muon.trackDist(i - 4, 2, arbitrationType);
       if (thisTrackDist < 999999) {  //current "raw" info that a track is close to a chamber
         ++nr_of_stations_crossed;
         station_was_crossed[i - 1] = 1;
         if (thisTrackDist > -10.)
           stations_w_track_at_boundary[i - 1] = thisTrackDist;
         else
           stations_w_track_at_boundary[i - 1] = 0.;
       }
       //current "raw" info that a segment is matched to the current track
       if (muon.segmentX(i - 4, 2, arbitrationType) < 999999) {
         station_has_segmentmatch[i - 1] = 1;
       }
     }
     // rough estimation of chamber border efficiency (should be parametrized better, this is just a quick guess):
     // TF1 * merf = new TF1("merf","-0.5*(TMath::Erf(x/6.)-1)",-100,100);
     // use above value to "unpunish" missing segment if close to border, i.e. rather than not adding any weight, add
     // the one from the function. Only for dist ~> -10 cm, else full punish!.
 
     // ********************************************************;
     // *** fill local info for this muon (do some counting) ***;
     // ************** end *************************************;
   }
 
   // ********************************************************;
   // *** calculate weights for each station *****************;
   // ************** begin ***********************************;
   //    const float slope = 0.5;
   //    const float attenuate_weight_regain = 1.;
   // if attenuate_weight_regain < 1., additional punishment if track is close to boundary and no segment
   const float attenuate_weight_regain = 0.5;
 
   for (int i = 1; i <= 8; ++i) {  // loop over all possible stations
 
     // first set all weights if a station has been crossed
     // later penalize if a station did not have a matching segment
 
     //old logic      if(station_has_segmentmatch[i-1] > 0 ) { // the track has an associated segment at the current station
     if (station_was_crossed[i - 1] > 0) {  // the track crossed this chamber (or was nearby)
       // - Apply a weight depending on the "depth" of the muon passage.
       // - The station_weight is later reduced for stations with badly matched segments.
       // - Even if there is no segment but the track passes close to a chamber boundary, the
       //   weight is set non zero and can go up to 0.5 of the full weight if the track is quite
       //   far from any station.
       ++position_in_stations;
 
       switch (nr_of_stations_crossed) {  // define different weights depending on how many stations were crossed
         case 1:
           station_weight[i - 1] = 1.f;
           break;
         case 2:
           if (position_in_stations == 1)
             station_weight[i - 1] = 0.33f;
           else
             station_weight[i - 1] = 0.67f;
           break;
         case 3:
           if (position_in_stations == 1)
             station_weight[i - 1] = 0.23f;
           else if (position_in_stations == 2)
             station_weight[i - 1] = 0.33f;
           else
             station_weight[i - 1] = 0.44f;
           break;
         case 4:
           if (position_in_stations == 1)
             station_weight[i - 1] = 0.10f;
           else if (position_in_stations == 2)
             station_weight[i - 1] = 0.20f;
           else if (position_in_stations == 3)
             station_weight[i - 1] = 0.30f;
           else
             station_weight[i - 1] = 0.40f;
           break;
 
         default:
           //    LogTrace("MuonIdentification")<<"            // Message: A muon candidate track has more than 4 stations with matching segments.";
           //    LogTrace("MuonIdentification")<<"            // Did not expect this - please let me know: ibloch@fnal.gov";
           // for all other cases
           station_weight[i - 1] = 1.f / nr_of_stations_crossed;
       }
 
       if (use_weight_regain_at_chamber_boundary) {  // reconstitute some weight if there is no match but the segment is close to a boundary:
         if (station_has_segmentmatch[i - 1] <= 0 && stations_w_track_at_boundary[i - 1] != 0.) {
           // if segment is not present but track in inefficient region, do not count as "missing match" but add some reduced weight.
           // original "match weight" is currently reduced by at least attenuate_weight_regain, variing with an error function down to 0 if the track is
           // inside the chamber.
           // remark: the additional scale of 0.5 normalizes Err to run from 0 to 1 in y
           station_weight[i - 1] = station_weight[i - 1] * attenuate_weight_regain * 0.5f *
                                   (std::erf(stations_w_track_at_boundary[i - 1] / 6.f) + 1.f);
         } else if (station_has_segmentmatch[i - 1] <= 0 &&
                    stations_w_track_at_boundary[i - 1] == 0.f) {  // no segment match and track well inside chamber
           // full penalization
           station_weight[i - 1] = 0.f;
         }
       } else {  // always fully penalize tracks with no matching segment, whether the segment is close to the boundary or not.
         if (station_has_segmentmatch[i - 1] <= 0)
           station_weight[i - 1] = 0.f;
       }
 
       // if track has matching segment, but the matching is not high quality, penalize
       if (station_has_segmentmatch[i - 1] > 0) {
         if (i <= 4) {  // we are in the DTs
           if (muon.dY(i, 1, arbitrationType) < 999999.f &&
               muon.dX(i, 1, arbitrationType) < 999999.f) {  // have both X and Y match
             const float pullTot2 =
                 std::pow(muon.pullX(i, 1, arbitrationType), 2.) + std::pow(muon.pullY(i, 1, arbitrationType), 2.);
             if (pullTot2 > 1.f) {
               const float dxy2 =
                   std::pow(muon.dX(i, 1, arbitrationType), 2.) + std::pow(muon.dY(i, 1, arbitrationType), 2.);
               // reduce weight
               if (use_match_dist_penalty) {
                 // only use pull if 3 sigma is not smaller than 3 cm
                 if (dxy2 < 9.f && pullTot2 > 9.f) {
                   if (dxy2 > 1.f)
                     station_weight[i - 1] *= 1.f / std::pow(dxy2, .125);
                 } else {
                   station_weight[i - 1] *= 1.f / std::pow(pullTot2, .125);
                 }
               }
             }
           } else if (muon.dY(i, 1, arbitrationType) >= 999999.f) {  // has no match in Y
             // has a match in X. Pull larger that 1 to avoid increasing the weight (just penalize, don't anti-penalize)
             if (muon.pullX(i, 1, arbitrationType) > 1.f) {
               // reduce weight
               if (use_match_dist_penalty) {
                 // only use pull if 3 sigma is not smaller than 3 cm
                 if (muon.dX(i, 1, arbitrationType) < 3.f && muon.pullX(i, 1, arbitrationType) > 3.f) {
                   if (muon.dX(i, 1, arbitrationType) > 1.f)
                     station_weight[i - 1] *= 1.f / std::pow(muon.dX(i, 1, arbitrationType), .25);
                 } else {
                   station_weight[i - 1] *= 1.f / std::pow(muon.pullX(i, 1, arbitrationType), .25);
                 }
               }
             }
           } else {  // has no match in X
             // has a match in Y. Pull larger that 1 to avoid increasing the weight (just penalize, don't anti-penalize)
             if (muon.pullY(i, 1, arbitrationType) > 1.f) {
               // reduce weight
               if (use_match_dist_penalty) {
                 // only use pull if 3 sigma is not smaller than 3 cm
                 if (muon.dY(i, 1, arbitrationType) < 3. && muon.pullY(i, 1, arbitrationType) > 3.) {
                   if (muon.dY(i, 1, arbitrationType) > 1.f)
                     station_weight[i - 1] *= 1.f / std::pow(muon.dY(i, 1, arbitrationType), .25);
                 } else {
                   station_weight[i - 1] *= 1.f / std::pow(muon.pullY(i, 1, arbitrationType), .25);
                 }
               }
             }
           }
         } else {  // We are in the CSCs
           const float pullTot2 =
               std::pow(muon.pullX(i - 4, 2, arbitrationType), 2.) + std::pow(muon.pullY(i - 4, 2, arbitrationType), 2.);
           if (pullTot2 > 1.f) {
             // reduce weight
             if (use_match_dist_penalty) {
               const float dxy2 =
                   std::pow(muon.dX(i - 4, 2, arbitrationType), 2.) + std::pow(muon.dY(i - 4, 2, arbitrationType), 2.);
               // only use pull if 3 sigma is not smaller than 3 cm
               if (dxy2 < 9.f && pullTot2 > 9.f) {
                 if (dxy2 > 1.f)
                   station_weight[i - 1] *= 1.f / std::pow(dxy2, .125);
               } else {
                 station_weight[i - 1] *= 1.f / std::pow(pullTot2, .125);
               }
             }
           }
         }
       }
 
       // Thoughts:
       // - should penalize if the segment has only x OR y info
       // - should also use the segment direction, as it now works!
 
     } else {  // track did not pass a chamber in this station - just reset weight
       station_weight[i - 1] = 0.f;
     }
 
     //increment final weight for muon:
     full_weight += station_weight[i - 1];
   }
 
   // if we don't expect any matches, we set the compatibility to
   // 0.5 as the track is as compatible with a muon as it is with
   // background - we should maybe rather set it to -0.5!
   if (nr_of_stations_crossed == 0) {
     //      full_weight = attenuate_weight_regain*0.5;
     full_weight = 0.5f;
   }
 
   // ********************************************************;
   // *** calculate weights for each station *****************;
   // ************** end *************************************;
 
   return full_weight;
 }
 
 bool muon::isGoodMuon(const reco::Muon& muon,
                       AlgorithmType type,
                       double minCompatibility,
                       reco::Muon::ArbitrationType arbitrationType) {
   if (!muon.isMatchesValid())
     return false;
   bool goodMuon = false;
 
   switch (type) {
     case TM2DCompatibility:
       // Simplistic first cut in the 2D segment- vs calo-compatibility plane. Will have to be refined!
       if (((0.8 * caloCompatibility(muon)) + (1.2 * segmentCompatibility(muon, arbitrationType))) > minCompatibility)
         goodMuon = true;
       else
         goodMuon = false;
       return goodMuon;
       break;
     default:
       //    LogTrace("MuonIdentification")<<"            // Invalid Algorithm Type called!";
       goodMuon = false;
       return goodMuon;
       break;
   }
 }
 
 bool muon::isGoodMuon(const reco::Muon& muon,
                       AlgorithmType type,
                       int minNumberOfMatches,
                       double maxAbsDx,
                       double maxAbsPullX,
                       double maxAbsDy,
                       double maxAbsPullY,
                       double maxChamberDist,
                       double maxChamberDistPull,
                       reco::Muon::ArbitrationType arbitrationType,
                       bool syncMinNMatchesNRequiredStationsInBarrelOnly,
                       bool applyAlsoAngularCuts) {
   if (!muon.isMatchesValid())
     return false;
   bool goodMuon = false;
 
   if (type == TMLastStation) {
     // To satisfy my own paranoia, if the user specifies that the
     // minimum number of matches is zero, then return true.
     if (minNumberOfMatches == 0)
       return true;
 
     unsigned int theStationMask = muon.stationMask(arbitrationType);
     unsigned int theRequiredStationMask =
         RequiredStationMask(muon, maxChamberDist, maxChamberDistPull, arbitrationType);
 
     // Require that there be at least a minimum number of segments
     int numSegs = 0;
     int numRequiredStations = 0;
     for (int it = 0; it < 8; ++it) {
       if (theStationMask & 1 << it)
         ++numSegs;
       if (theRequiredStationMask & 1 << it)
         ++numRequiredStations;
     }
 
     // Make sure the minimum number of matches is not greater than
     // the number of required stations but still greater than zero
     if (syncMinNMatchesNRequiredStationsInBarrelOnly) {
       // Note that we only do this in the barrel region!
       if (std::abs(muon.eta()) < 1.2) {
         if (minNumberOfMatches > numRequiredStations)
           minNumberOfMatches = numRequiredStations;
         if (minNumberOfMatches < 1)  //SK: this only happens for negative values
           minNumberOfMatches = 1;
       }
     } else {
       if (minNumberOfMatches > numRequiredStations)
         minNumberOfMatches = numRequiredStations;
       if (minNumberOfMatches < 1)  //SK: this only happens for negative values
         minNumberOfMatches = 1;
     }
 
     if (numSegs >= minNumberOfMatches)
       goodMuon = true;
 
     // Require that last required station have segment
     // If there are zero required stations keep track
     // of the last station with a segment so that we may
     // apply the quality cuts below to it instead
     int lastSegBit = 0;
     if (theRequiredStationMask) {
       for (int stationIdx = 7; stationIdx >= 0; --stationIdx)
         if (theRequiredStationMask & 1 << stationIdx) {
           if (theStationMask & 1 << stationIdx) {
             lastSegBit = stationIdx;
             goodMuon &= 1;
             break;
           } else {
             goodMuon = false;
             break;
           }
         }
     } else {
       for (int stationIdx = 7; stationIdx >= 0; --stationIdx)
         if (theStationMask & 1 << stationIdx) {
           lastSegBit = stationIdx;
           break;
         }
     }
 
     if (!goodMuon)
       return false;
 
     // Impose pull cuts on last segment
     int station = 0, detector = 0;
     station = lastSegBit < 4 ? lastSegBit + 1 : lastSegBit - 3;
     detector = lastSegBit < 4 ? 1 : 2;
 
     // Check x information
     if (std::abs(muon.pullX(station, detector, arbitrationType, true)) > maxAbsPullX &&
         std::abs(muon.dX(station, detector, arbitrationType)) > maxAbsDx)
       return false;
 
     if (applyAlsoAngularCuts && std::abs(muon.pullDxDz(station, detector, arbitrationType, true)) > maxAbsPullX)
       return false;
 
     // Is this a tight algorithm, i.e. do we bother to check y information?
     if (maxAbsDy < 999999) {  // really if maxAbsDy < 1E9 as currently defined
 
       // Check y information
       if (detector == 2) {  // CSC
         if (std::abs(muon.pullY(station, 2, arbitrationType, true)) > maxAbsPullY &&
             std::abs(muon.dY(station, 2, arbitrationType)) > maxAbsDy)
           return false;
 
         if (applyAlsoAngularCuts && std::abs(muon.pullDyDz(station, 2, arbitrationType, true)) > maxAbsPullY)
           return false;
       } else {
         //
         // In DT, if this is a "Tight" algorithm and the last segment is
         // missing y information (always the case in station 4!!!), impose
         // respective cuts on the next station in the stationMask that has
         // a segment with y information.  If there are no segments with y
         // information then there is nothing to penalize. Should we
         // penalize in Tight for having zero segments with y information?
         // That is the fundamental question.  Of course I am being uber
         // paranoid; if this is a good muon then there will probably be at
         // least one segment with y information but not always.  Suppose
         // somehow a muon only creates segments in station 4, then we
         // definitely do not want to require that there be at least one
         // segment with y information because we will lose it completely.
         //
 
         for (int stationIdx = station; stationIdx > 0; --stationIdx) {
           if (!(theStationMask & 1 << (stationIdx - 1)))  // don't bother if the station is not in the stationMask
             continue;
 
           if (muon.dY(stationIdx, 1, arbitrationType) > 999998)  // no y-information
             continue;
 
           if (std::abs(muon.pullY(stationIdx, 1, arbitrationType, true)) > maxAbsPullY &&
               std::abs(muon.dY(stationIdx, 1, arbitrationType)) > maxAbsDy) {
             return false;
           }
 
           if (applyAlsoAngularCuts && std::abs(muon.pullDyDz(stationIdx, 1, arbitrationType, true)) > maxAbsPullY)
             return false;
 
           // If we get this far then great this is a good muon
           return true;
         }
       }
     }
 
     return goodMuon;
   }  // TMLastStation
 
   // TMOneStation requires only that there be one "good" segment, regardless
   // of the required stations.  We do not penalize if there are absolutely zero
   // segments with y information in the Tight algorithm.  Maybe I'm being
   // paranoid but so be it.  If it's really a good muon then we will probably
   // find at least one segment with both x and y information but you never
   // know, and I don't want to deal with a potential inefficiency in the DT
   // like we did with the original TMLastStation.  Incidentally, not penalizing
   // for total lack of y information in the Tight algorithm is what is done in
   // the new TMLastStation
   //
   if (type == TMOneStation) {
     unsigned int theStationMask = muon.stationMask(arbitrationType);
 
     // Of course there must be at least one segment
     if (!theStationMask)
       return false;
 
     int station = 0, detector = 0;
     // Keep track of whether or not there is a DT segment with y information.
     // In the end, if it turns out there are absolutely zero DT segments with
     // y information, require only that there was a segment with good x info.
     // This of course only applies to the Tight algorithms.
     bool existsGoodDTSegX = false;
     bool existsDTSegY = false;
 
     // Impose cuts on the segments in the station mask until we find a good one
     // Might as well start with the lowest bit to speed things up.
     for (int stationIdx = 0; stationIdx <= 7; ++stationIdx)
       if (theStationMask & 1 << stationIdx) {
         station = stationIdx < 4 ? stationIdx + 1 : stationIdx - 3;
         detector = stationIdx < 4 ? 1 : 2;
 
         if ((std::abs(muon.pullX(station, detector, arbitrationType, true)) > maxAbsPullX &&
              std::abs(muon.dX(station, detector, arbitrationType)) > maxAbsDx) ||
             (applyAlsoAngularCuts && std::abs(muon.pullDxDz(station, detector, arbitrationType, true)) > maxAbsPullX))
           continue;
         else if (detector == 1)
           existsGoodDTSegX = true;
 
         // Is this a tight algorithm?  If yes, use y information
         if (maxAbsDy < 999999) {
           if (detector == 2) {  // CSC
             if ((std::abs(muon.pullY(station, 2, arbitrationType, true)) > maxAbsPullY &&
                  std::abs(muon.dY(station, 2, arbitrationType)) > maxAbsDy) ||
                 (applyAlsoAngularCuts && std::abs(muon.pullDyDz(station, 2, arbitrationType, true)) > maxAbsPullY))
               continue;
           } else {
             if (muon.dY(station, 1, arbitrationType) > 999998)  // no y-information
               continue;
             else
               existsDTSegY = true;
 
             if ((std::abs(muon.pullY(station, 1, arbitrationType, true)) > maxAbsPullY &&
                  std::abs(muon.dY(station, 1, arbitrationType)) > maxAbsDy) ||
                 (applyAlsoAngularCuts && std::abs(muon.pullDyDz(station, 1, arbitrationType, true)) > maxAbsPullY)) {
               continue;
             }
           }
         }
 
         // If we get this far then great this is a good muon
         return true;
       }
 
     // If we get this far then for sure there are no "good" CSC segments. For
     // DT, check if there were any segments with y information.  If there
     // were none, but there was a segment with good x, then we're happy. If
     // there WERE segments with y information, then they must have been shit
     // since we are here so fail it.  Of course, if this is a Loose algorithm
     // then fail immediately since if we had good x we would already have
     // returned true
     if (maxAbsDy < 999999) {
       if (existsDTSegY)
         return false;
       else if (existsGoodDTSegX)
         return true;
     } else
       return false;
   }  // TMOneStation
 
   if (type == RPCMu) {
     if (minNumberOfMatches == 0)
       return true;
 
     int nMatch = 0;
     for (const auto& chamberMatch : muon.matches()) {
       if (chamberMatch.detector() != MuonSubdetId::RPC)
         continue;
 
       const float trkX = chamberMatch.x;
       const float errX = chamberMatch.xErr;
 
       for (const auto& rpcMatch : chamberMatch.rpcMatches) {
         const float rpcX = rpcMatch.x;
         const float dX = std::abs(rpcX - trkX);
         if (dX < maxAbsDx or dX < maxAbsPullX * errX) {
           ++nMatch;
           break;
         }
       }
     }
 
     if (nMatch >= minNumberOfMatches)
       return true;
     else
       return false;
   }  // RPCMu
 
   if (type == ME0Mu) {
     if (minNumberOfMatches == 0)
       return true;
 
     int nMatch = 0;
     for (const auto& chamberMatch : muon.matches()) {
       if (chamberMatch.detector() != MuonSubdetId::ME0)
         continue;
 
       const float trkX = chamberMatch.x;
       const float errX2 = chamberMatch.xErr * chamberMatch.xErr;
       const float trkY = chamberMatch.y;
       const float errY2 = chamberMatch.yErr * chamberMatch.yErr;
 
       for (const auto& segment : chamberMatch.me0Matches) {
         const float me0X = segment.x;
         const float me0ErrX2 = segment.xErr * segment.xErr;
         const float me0Y = segment.y;
         const float me0ErrY2 = segment.yErr * segment.yErr;
 
         const float dX = std::abs(me0X - trkX);
         const float dY = std::abs(me0Y - trkY);
         const float invPullX2 = errX2 + me0ErrX2;
         const float invPullY2 = errY2 + me0ErrY2;
 
         if ((dX < maxAbsDx or dX < maxAbsPullX * std::sqrt(invPullX2)) and
             (dY < maxAbsDy or dY < maxAbsPullY * std::sqrt(invPullY2))) {
           ++nMatch;
           break;
         }
       }
     }
 
     return (nMatch >= minNumberOfMatches);
   }  // ME0Mu
 
   if (type == GEMMu) {
     if (minNumberOfMatches == 0)
       return true;
 
     int nMatch = 0;
     for (const auto& chamberMatch : muon.matches()) {
       if (chamberMatch.detector() != MuonSubdetId::GEM)
         continue;
 
       const float trkX = chamberMatch.x;
       const float errX2 = chamberMatch.xErr * chamberMatch.xErr;
       const float trkY = chamberMatch.y;
       const float errY2 = chamberMatch.yErr * chamberMatch.yErr;
 
       for (const auto& segment : chamberMatch.gemMatches) {
         const float gemX = segment.x;
         const float gemErrX2 = segment.xErr * segment.xErr;
         const float gemY = segment.y;
         const float gemErrY2 = segment.yErr * segment.yErr;
 
         const float dX = std::abs(gemX - trkX);
         const float dY = std::abs(gemY - trkY);
         const float invPullX2 = errX2 + gemErrX2;
         const float invPullY2 = errY2 + gemErrY2;
 
         if ((dX < maxAbsDx or dX < maxAbsPullX * std::sqrt(invPullX2)) and
             (dY < maxAbsDy or dY < maxAbsPullY * std::sqrt(invPullY2))) {
           ++nMatch;
           break;
         }
       }
     }
 
     return (nMatch >= minNumberOfMatches);
   }  // GEMMu
 
   return goodMuon;
 }
 
 bool muon::isGoodMuon(const reco::Muon& muon, SelectionType type, reco::Muon::ArbitrationType arbitrationType) {
   switch (type) {
     case muon::All:
       return true;
       break;
     case muon::AllGlobalMuons:
       return muon.isGlobalMuon();
       break;
     case muon::AllTrackerMuons:
       return muon.isTrackerMuon();
       break;
     case muon::AllStandAloneMuons:
       return muon.isStandAloneMuon();
       break;
     case muon::TrackerMuonArbitrated:
       return muon.isTrackerMuon() && muon.numberOfMatches(arbitrationType) > 0;
       break;
     case muon::AllArbitrated:
       return !muon.isTrackerMuon() || muon.numberOfMatches(arbitrationType) > 0;
       break;
     case muon::GlobalMuonPromptTight:
       return muon.isGlobalMuon() && muon.globalTrack()->normalizedChi2() < 10. &&
              muon.globalTrack()->hitPattern().numberOfValidMuonHits() > 0;
       break;
       // For "Loose" algorithms we choose maximum y quantity cuts of 1E9 instead of
       // 9999 as before.  We do this because the muon methods return 999999 (note
       // there are six 9's) when the requested information is not available.  For
       // example, if a muon fails to traverse the z measuring superlayer in a station
       // in the DT, then all methods involving segmentY in this station return
       // 999999 to demonstrate that the information is missing.  In order to not
       // penalize muons for missing y information in Loose algorithms where we do
       // not care at all about y information, we raise these limits.  In the
       // TMLastStation and TMOneStation algorithms we actually use this huge number
       // to determine whether to consider y information at all.
     case muon::TMLastStationLoose:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, true, false);
       break;
     case muon::TMLastStationTight:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, true, false);
       break;
     case muon::TMOneStationLoose:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, false);
       break;
     case muon::TMOneStationTight:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, false);
       break;
     case muon::TMLastStationOptimizedLowPtLoose:
       if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, false);
       else
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, false, false);
       break;
     case muon::TMLastStationOptimizedLowPtTight:
       if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, false);
       else
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, false, false);
       break;
       //compatibility loose
     case muon::TM2DCompatibilityLoose:
       return muon.isTrackerMuon() && isGoodMuon(muon, TM2DCompatibility, 0.7, arbitrationType);
       break;
       //compatibility tight
     case muon::TM2DCompatibilityTight:
       return muon.isTrackerMuon() && isGoodMuon(muon, TM2DCompatibility, 1.0, arbitrationType);
       break;
     case muon::GMTkChiCompatibility:
       return muon.isGlobalMuon() && muon.isQualityValid() &&
              std::abs(muon.combinedQuality().trkRelChi2 - muon.innerTrack()->normalizedChi2()) < 2.0;
       break;
     case muon::GMStaChiCompatibility:
       return muon.isGlobalMuon() && muon.isQualityValid() &&
              std::abs(muon.combinedQuality().staRelChi2 - muon.outerTrack()->normalizedChi2()) < 2.0;
       break;
     case muon::GMTkKinkTight:
       return muon.isGlobalMuon() && muon.isQualityValid() && muon.combinedQuality().trkKink < 100.0;
       break;
     case muon::TMLastStationAngLoose:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, false, true);
       break;
     case muon::TMLastStationAngTight:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, false, true);
       break;
     case muon::TMOneStationAngLoose:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, true);
       break;
     case muon::TMOneStationAngTight:
       return muon.isTrackerMuon() &&
              isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, true);
       break;
     case muon::TMLastStationOptimizedBarrelLowPtLoose:
       if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMOneStation, 1, 3, 3, 1E9, 1E9, 1E9, 1E9, arbitrationType, false, false);
       else
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMLastStation, 2, 3, 3, 1E9, 1E9, -3, -3, arbitrationType, true, false);
       break;
     case muon::TMLastStationOptimizedBarrelLowPtTight:
       if (muon.pt() < 8. && std::abs(muon.eta()) < 1.2)
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMOneStation, 1, 3, 3, 3, 3, 1E9, 1E9, arbitrationType, false, false);
       else
         return muon.isTrackerMuon() &&
                isGoodMuon(muon, TMLastStation, 2, 3, 3, 3, 3, -3, -3, arbitrationType, true, false);
       break;
     case muon::RPCMuLoose:
       return muon.isRPCMuon() && isGoodMuon(muon, RPCMu, 2, 20, 4, 1e9, 1e9, 1e9, 1e9, arbitrationType, false, false);
       break;
     case muon::AllME0Muons:
       return muon.isME0Muon();
       break;
     case muon::ME0MuonArbitrated:
       return muon.isME0Muon() &&
              isGoodMuon(muon, ME0Mu, 1, 1e9, 1e9, 1e9, 1e9, 1e9, 1e9, arbitrationType, false, false);
       break;
     case muon::AllGEMMuons:
       return muon.isGEMMuon();
       break;
     case muon::GEMMuonArbitrated:
       return muon.isGEMMuon() &&
              isGoodMuon(muon, GEMMu, 1, 1e9, 1e9, 1e9, 1e9, 1e9, 1e9, arbitrationType, false, false);
       break;
     case muon::TriggerIdLoose:
       return isLooseTriggerMuon(muon);
       break;
     default:
       return false;
   }
 }
 
 bool muon::overlap(
     const reco::Muon& muon1, const reco::Muon& muon2, double pullX, double pullY, bool checkAdjacentChambers) {
   unsigned int nMatches1 = muon1.numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
   unsigned int nMatches2 = muon2.numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
   unsigned int betterMuon = (muon1.pt() > muon2.pt() ? 1 : 2);
   for (std::vector<reco::MuonChamberMatch>::const_iterator chamber1 = muon1.matches().begin();
        chamber1 != muon1.matches().end();
        ++chamber1)
     for (std::vector<reco::MuonChamberMatch>::const_iterator chamber2 = muon2.matches().begin();
          chamber2 != muon2.matches().end();
          ++chamber2) {
       // if ( (chamber1->segmentMatches.empty() || chamber2->segmentMatches.empty()) ) continue;
 
       // handle case where both muons have information about the same chamber
       // here we know how close they are
       if (chamber1->id == chamber2->id) {
         // found the same chamber
         if (std::abs(chamber1->x - chamber2->x) <
             pullX * sqrt(chamber1->xErr * chamber1->xErr + chamber2->xErr * chamber2->xErr)) {
           if (betterMuon == 1)
             nMatches2--;
           else
             nMatches1--;
           if (nMatches1 == 0 || nMatches2 == 0)
             return true;
           continue;
         }
         if (std::abs(chamber1->y - chamber2->y) <
             pullY * sqrt(chamber1->yErr * chamber1->yErr + chamber2->yErr * chamber2->yErr)) {
           if (betterMuon == 1)
             nMatches2--;
           else
             nMatches1--;
           if (nMatches1 == 0 || nMatches2 == 0)
             return true;
         }
       } else {
         if (!checkAdjacentChambers)
           continue;
         // check if tracks are pointing into overlaping region of the CSC detector
         if (chamber1->id.subdetId() != MuonSubdetId::CSC || chamber2->id.subdetId() != MuonSubdetId::CSC)
           continue;
         CSCDetId id1(chamber1->id);
         CSCDetId id2(chamber2->id);
         if (id1.endcap() != id2.endcap())
           continue;
         if (id1.station() != id2.station())
           continue;
         if (id1.ring() != id2.ring())
           continue;
         if (std::abs(id1.chamber() - id2.chamber()) > 1)
           continue;
         // FIXME: we don't handle 18->1; 36->1 transitions since
         // I don't know how to check for sure how many chambers
         // are there. Probably need to hard code some checks.
 
         // Now we have to make sure that both tracks are close to an edge
         // FIXME: ignored Y coordinate for now
         if (std::abs(chamber1->edgeX) > chamber1->xErr * pullX)
           continue;
         if (std::abs(chamber2->edgeX) > chamber2->xErr * pullX)
           continue;
         if (chamber1->x * chamber2->x < 0) {  // check if the same edge
           if (betterMuon == 1)
             nMatches2--;
           else
             nMatches1--;
           if (nMatches1 == 0 || nMatches2 == 0)
             return true;
         }
       }
     }
   return false;
 }
 
 bool muon::isLooseTriggerMuon(const reco::Muon& muon) {
   // Requirements:
   // - no depencence on information not availabe in the muon object
   // - use only robust inputs
   bool tk_id = muon::isGoodMuon(muon, TMOneStationTight);
   if (not tk_id)
     return false;
   bool layer_requirements = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
                             muon.innerTrack()->hitPattern().pixelLayersWithMeasurement() > 0;
   bool match_requirements =
       (muon.expectedNnumberOfMatchedStations() < 2) or (muon.numberOfMatchedStations() > 1) or (muon.pt() < 8);
   return layer_requirements and match_requirements;
 }
 
 bool muon::isTightMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
   if (!muon.isPFMuon() || !muon.isGlobalMuon())
     return false;
 
   bool muID = isGoodMuon(muon, GlobalMuonPromptTight) && (muon.numberOfMatchedStations() > 1);
 
   bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
               muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0;
 
   bool ip = std::abs(muon.muonBestTrack()->dxy(vtx.position())) < 0.2 &&
             std::abs(muon.muonBestTrack()->dz(vtx.position())) < 0.5;
 
   return muID && hits && ip;
 }
 
 bool muon::isLooseMuon(const reco::Muon& muon) {
   return muon.isPFMuon() && (muon.isGlobalMuon() || muon.isTrackerMuon());
 }
 
 bool muon::isMediumMuon(const reco::Muon& muon, bool run2016_hip_mitigation) {
   if (not isLooseMuon(muon))
     return false;
   if (run2016_hip_mitigation) {
     if (muon.innerTrack()->validFraction() < 0.49)
       return false;
   } else {
     if (muon.innerTrack()->validFraction() < 0.8)
       return false;
   }
 
   bool goodGlb = muon.isGlobalMuon() && muon.globalTrack()->normalizedChi2() < 3. &&
                  muon.combinedQuality().chi2LocalPosition < 12. && muon.combinedQuality().trkKink < 20.;
 
   return (segmentCompatibility(muon) > (goodGlb ? 0.303 : 0.451));
 }
 
 bool muon::isSoftMuon(const reco::Muon& muon, const reco::Vertex& vtx, bool run2016_hip_mitigation) {
   bool muID = muon::isGoodMuon(muon, TMOneStationTight);
 
   if (!muID)
     return false;
 
   bool layers = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
                 muon.innerTrack()->hitPattern().pixelLayersWithMeasurement() > 0;
 
   bool ishighq = muon.innerTrack()->quality(reco::Track::highPurity);
 
   bool ip =
       std::abs(muon.innerTrack()->dxy(vtx.position())) < 0.3 && std::abs(muon.innerTrack()->dz(vtx.position())) < 20.;
 
   return layers && ip && (ishighq | run2016_hip_mitigation);
 }
 
 bool muon::isHighPtMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
   if (!muon.isGlobalMuon())
     return false;
 
   bool muValHits = (muon.globalTrack()->hitPattern().numberOfValidMuonHits() > 0 ||
                     muon.tunePMuonBestTrack()->hitPattern().numberOfValidMuonHits() > 0);
 
   bool muMatchedSt = muon.numberOfMatchedStations() > 1;
   if (!muMatchedSt) {
     if (muon.isTrackerMuon() && muon.numberOfMatchedStations() == 1) {
       if (muon.expectedNnumberOfMatchedStations() < 2 || !(muon.stationMask() == 1 || muon.stationMask() == 16) ||
           muon.numberOfMatchedRPCLayers() > 2)
         muMatchedSt = true;
     }
   }
 
   bool muID = muValHits && muMatchedSt;
 
   bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
               muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0;
 
   bool momQuality = muon.tunePMuonBestTrack()->ptError() / muon.tunePMuonBestTrack()->pt() < 0.3;
 
   bool ip =
       std::abs(muon.innerTrack()->dxy(vtx.position())) < 0.2 && std::abs(muon.innerTrack()->dz(vtx.position())) < 0.5;
 
   return muID && hits && momQuality && ip;
 }
 
 bool muon::isTrackerHighPtMuon(const reco::Muon& muon, const reco::Vertex& vtx) {
   bool muID = muon.isTrackerMuon() && muon.track().isNonnull() && (muon.numberOfMatchedStations() > 1);
   if (!muID)
     return false;
 
   bool hits = muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
               muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0;
 
   bool momQuality = muon.tunePMuonBestTrack()->ptError() < 0.3 * muon.tunePMuonBestTrack()->pt();
 
   bool ip =
       std::abs(muon.innerTrack()->dxy(vtx.position())) < 0.2 && std::abs(muon.innerTrack()->dz(vtx.position())) < 0.5;
 
   return muID && hits && momQuality && ip;
 }
 
 int muon::sharedSegments(const reco::Muon& mu, const reco::Muon& mu2, unsigned int segmentArbitrationMask) {
   int ret = 0;
 
   // Will do with a stupid double loop, since creating and filling a map is probably _more_ inefficient for a single lookup.
   for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = mu.matches().begin();
        chamberMatch != mu.matches().end();
        ++chamberMatch) {
     if (chamberMatch->segmentMatches.empty())
       continue;
     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch2 = mu2.matches().begin();
          chamberMatch2 != mu2.matches().end();
          ++chamberMatch2) {
       if (chamberMatch2->segmentMatches.empty())
         continue;
       if (chamberMatch2->id() != chamberMatch->id())
         continue;
       for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();
            segmentMatch != chamberMatch->segmentMatches.end();
            ++segmentMatch) {
         if (!segmentMatch->isMask(segmentArbitrationMask))
           continue;
         for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch2 = chamberMatch2->segmentMatches.begin();
              segmentMatch2 != chamberMatch2->segmentMatches.end();
              ++segmentMatch2) {
           if (!segmentMatch2->isMask(segmentArbitrationMask))
             continue;
           if ((segmentMatch->cscSegmentRef.isNonnull() &&
                segmentMatch->cscSegmentRef == segmentMatch2->cscSegmentRef) ||
               (segmentMatch->dtSegmentRef.isNonnull() && segmentMatch->dtSegmentRef == segmentMatch2->dtSegmentRef)) {
             ++ret;
           }  // is the same
         }    // segment of mu2 in chamber
       }      // segment of mu1 in chamber
     }        // chamber of mu2
   }          // chamber of mu1
 
   return ret;
 }
 
 bool outOfTimeMuon(const reco::Muon& muon) {
   const auto& combinedTime = muon.time();
   const auto& rpcTime = muon.rpcTime();
   bool combinedTimeIsOk = (combinedTime.nDof > 7);
   bool rpcTimeIsOk = (rpcTime.nDof > 1 && std::abs(rpcTime.timeAtIpInOutErr) < 0.001);
   bool outOfTime = false;
   if (rpcTimeIsOk) {
     if ((std::abs(rpcTime.timeAtIpInOut) > 10) && !(combinedTimeIsOk && std::abs(combinedTime.timeAtIpInOut) < 10))
       outOfTime = true;
   } else {
     if (combinedTimeIsOk && (combinedTime.timeAtIpInOut > 20 || combinedTime.timeAtIpInOut < -45))
       outOfTime = true;
   }
   return outOfTime;
 }
 
 reco::Muon::Selector muon::makeSelectorBitset(reco::Muon const& muon,
                                               reco::Vertex const* vertex,
                                               bool run2016_hip_mitigation) {
   // https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideMuonIdRun2
   unsigned int selectors = muon.selectors();
   // Compute Id and Isolation variables
   double chIso = muon.pfIsolationR04().sumChargedHadronPt;
   double nIso = muon.pfIsolationR04().sumNeutralHadronEt;
   double phoIso = muon.pfIsolationR04().sumPhotonEt;
   double puIso = muon.pfIsolationR04().sumPUPt;
   double dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5 * puIso, 0.);
   double dbCorrectedRelIso = dbCorrectedIsolation / muon.pt();
   double tkRelIso = muon.isolationR03().sumPt / muon.pt();
 
   // Base selectors
   if (muon::isLooseMuon(muon))
     selectors |= (1UL << reco::Muon::CutBasedIdLoose);
   if (vertex) {
     if (muon::isTightMuon(muon, *vertex))
       selectors |= (1UL << reco::Muon::CutBasedIdTight);
     if (muon::isSoftMuon(muon, *vertex, run2016_hip_mitigation))
       selectors |= (1UL << reco::Muon::SoftCutBasedId);
     if (muon::isHighPtMuon(muon, *vertex))
       selectors |= (1UL << reco::Muon::CutBasedIdGlobalHighPt);
     if (muon::isTrackerHighPtMuon(muon, *vertex))
       selectors |= (1UL << reco::Muon::CutBasedIdTrkHighPt);
   }
   if (muon::isMediumMuon(muon, run2016_hip_mitigation)) {
     selectors |= (1UL << reco::Muon::CutBasedIdMedium);
     if (vertex and std::abs(muon.muonBestTrack()->dz(vertex->position())) < 0.1 and
         std::abs(muon.muonBestTrack()->dxy(vertex->position())) < 0.02)
       selectors |= (1UL << reco::Muon::CutBasedIdMediumPrompt);
   }
 
   // PF isolation
   if (dbCorrectedRelIso < 0.40)
     selectors |= (1UL << reco::Muon::PFIsoVeryLoose);
   if (dbCorrectedRelIso < 0.25)
     selectors |= (1UL << reco::Muon::PFIsoLoose);
   if (dbCorrectedRelIso < 0.20)
     selectors |= (1UL << reco::Muon::PFIsoMedium);
   if (dbCorrectedRelIso < 0.15)
     selectors |= (1UL << reco::Muon::PFIsoTight);
   if (dbCorrectedRelIso < 0.10)
     selectors |= (1UL << reco::Muon::PFIsoVeryTight);
   if (dbCorrectedRelIso < 0.05)
     selectors |= (1UL << reco::Muon::PFIsoVeryVeryTight);
 
   // Tracker isolation
   if (tkRelIso < 0.10)
     selectors |= (1UL << reco::Muon::TkIsoLoose);
   if (tkRelIso < 0.05)
     selectors |= (1UL << reco::Muon::TkIsoTight);
 
   // Trigger selectors
   if (isLooseTriggerMuon(muon))
     selectors |= (1UL << reco::Muon::TriggerIdLoose);
 
   // Timing
   if (!outOfTimeMuon(muon))
     selectors |= (1UL << reco::Muon::InTimeMuon);
 
   return static_cast<reco::Muon::Selector>(selectors);
 }

This page was automatically generated by the 2.2.1 LXR engine. The LXR team