diff --git a/DataFormats/MuonReco/interface/Muon.h b/DataFormats/MuonReco/interface/Muon.h
index 2f22dc6b2f4ac..78309ecf475ab 100755
--- a/DataFormats/MuonReco/interface/Muon.h
+++ b/DataFormats/MuonReco/interface/Muon.h
@@ -30,7 +30,7 @@ namespace reco {
     /// constructor from values
     Muon(  Charge, const LorentzVector &, const Point & = Point( 0, 0, 0 ) );
     /// create a clone
-    Muon * clone() const;
+    Muon * clone() const override;
 
     
     
@@ -46,20 +46,20 @@ namespace reco {
     /// reference to Track reconstructed in the tracker only
     using reco::RecoCandidate::track;
     virtual TrackRef innerTrack() const { return innerTrack_; }
-    virtual TrackRef track() const { return innerTrack(); }
+    TrackRef track() const override { return innerTrack(); }
     /// reference to Track reconstructed in the muon detector only
     virtual TrackRef outerTrack() const { return outerTrack_; }
-    virtual TrackRef standAloneMuon() const { return outerTrack(); }
+    TrackRef standAloneMuon() const override { return outerTrack(); }
     /// reference to Track reconstructed in both tracked and muon detector
     virtual TrackRef globalTrack() const { return globalTrack_; }
-    virtual TrackRef combinedMuon() const { return globalTrack(); }
+    TrackRef combinedMuon() const override { return globalTrack(); }
 
     virtual TrackRef tpfmsTrack() const { return muonTrackFromMap(TPFMS);}
     virtual TrackRef pickyTrack() const { return muonTrackFromMap(Picky);}
     virtual TrackRef dytTrack()   const { return muonTrackFromMap(DYT);}
     
-    virtual const Track * bestTrack() const         {return muonTrack(bestTrackType_).get();}
-    virtual TrackBaseRef  bestTrackRef() const      {return reco::TrackBaseRef(muonTrack(bestTrackType_));}
+    const Track * bestTrack() const override         {return muonTrack(bestTrackType_).get();}
+    TrackBaseRef  bestTrackRef() const override      {return reco::TrackBaseRef(muonTrack(bestTrackType_));}
     virtual TrackRef      muonBestTrack() const     {return muonTrack(bestTrackType_);}
     virtual MuonTrackType muonBestTrackType() const {return bestTrackType_;}
     virtual TrackRef      tunePMuonBestTrack() const     {return muonTrack(bestTunePTrackType_);}
@@ -182,6 +182,45 @@ namespace reco {
     enum ArbitrationType { NoArbitration, SegmentArbitration, SegmentAndTrackArbitration, SegmentAndTrackArbitrationCleaned,
 			   RPCHitAndTrackArbitration, GEMSegmentAndTrackArbitration, ME0SegmentAndTrackArbitration };
     
+    ///
+    /// ====================== STANDARD SELECTORS ===========================
+    ///
+    enum Selector {
+      CutBasedIdLoose        = 1UL<< 0,  
+      CutBasedIdMedium       = 1UL<< 1,  
+      CutBasedIdMediumPrompt = 1UL<< 2,  // medium with IP cuts
+      CutBasedIdTight        = 1UL<< 3,  
+      CutBasedIdGlobalHighPt = 1UL<< 4,  // high pt muon for Z',W' (better momentum resolution)
+      CutBasedIdTrkHighPt    = 1UL<< 5,  // high pt muon for boosted Z (better efficiency)
+      PFIsoVeryLoose         = 1UL<< 6,  // reliso<0.40
+      PFIsoLoose             = 1UL<< 7,  // reliso<0.25
+      PFIsoMedium            = 1UL<< 8,  // reliso<0.20
+      PFIsoTight             = 1UL<< 9,  // reliso<0.15
+      PFIsoVeryTight         = 1UL<<10,  // reliso<0.10
+      TkIsoLoose             = 1UL<<11,  // reliso<0.10
+      TkIsoTight             = 1UL<<12,  // reliso<0.05
+      SoftCutBasedId         = 1UL<<13,  
+      SoftMvaId              = 1UL<<14,  
+      MvaLoose               = 1UL<<15,  
+      MvaMedium              = 1UL<<16,  
+      MvaTight               = 1UL<<17,
+      MiniIsoLoose           = 1UL<<18,  // reliso<0.40
+      MiniIsoMedium          = 1UL<<19,  // reliso<0.20
+      MiniIsoTight           = 1UL<<20,  // reliso<0.10
+      MiniIsoVeryTight       = 1UL<<21   // reliso<0.05
+
+    };
+    
+    bool passed( unsigned int selection ) const { return (selectors_ & selection)==selection; }
+    unsigned int selectors() const { return selectors_; }
+    void setSelectors( unsigned int selectors ){ selectors_ = selectors; }
+    void setSelector(Selector selector, bool passed){ 
+      if (passed)
+	selectors_ |= selector;
+      else
+	selectors_ &= ~selector;
+    }
+
     ///
     /// ====================== USEFUL METHODs ===========================
     ///
@@ -221,11 +260,11 @@ namespace reco {
     void setType( unsigned int type ) { type_ = type; }
     unsigned int type() const { return type_; }
     // override of method in base class reco::Candidate
-    bool isMuon() const { return true; }
-    bool isGlobalMuon()     const { return type_ & GlobalMuon; }
-    bool isTrackerMuon()    const { return type_ & TrackerMuon; }
-    bool isStandAloneMuon() const { return type_ & StandAloneMuon; }
-    bool isCaloMuon() const { return type_ & CaloMuon; }
+    bool isMuon() const override { return true; }
+    bool isGlobalMuon()     const override { return type_ & GlobalMuon; }
+    bool isTrackerMuon()    const override { return type_ & TrackerMuon; }
+    bool isStandAloneMuon() const override { return type_ & StandAloneMuon; }
+    bool isCaloMuon() const override { return type_ & CaloMuon; }
     bool isPFMuon() const {return type_ & PFMuon;} //fix me ! Has to go to type
     bool isRPCMuon() const {return type_ & RPCMuon;}
     bool isGEMMuon() const {return type_ & GEMMuon;}
@@ -233,7 +272,7 @@ namespace reco {
     
   private:
     /// check overlap with another candidate
-    virtual bool overlap( const Candidate & ) const;
+    bool overlap( const Candidate & ) const override;
     /// reference to Track reconstructed in the tracker only
     TrackRef innerTrack_;
     /// reference to Track reconstructed in the muon detector only
@@ -288,7 +327,8 @@ namespace reco {
     /// get pointers to best segment and corresponding chamber in vector of chambers
     std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> pair( const std::vector<const MuonChamberMatch*> &,
 									ArbitrationType type = SegmentAndTrackArbitration ) const;
-     
+    /// selector bitmap
+    unsigned int selectors_;
    public:
      /// get number of segments
      int numberOfSegments( int station, int muonSubdetId, ArbitrationType type = SegmentAndTrackArbitration ) const;
diff --git a/DataFormats/MuonReco/interface/MuonSelectors.h b/DataFormats/MuonReco/interface/MuonSelectors.h
index df2501c026ccd..3a496015db9a9 100644
--- a/DataFormats/MuonReco/interface/MuonSelectors.h
+++ b/DataFormats/MuonReco/interface/MuonSelectors.h
@@ -87,9 +87,10 @@ namespace muon {
 
    bool isTightMuon(const reco::Muon&, const reco::Vertex&);
    bool isLooseMuon(const reco::Muon&);
-   bool isMediumMuon(const reco::Muon&);
-   bool isSoftMuon(const reco::Muon&, const reco::Vertex&);
+   bool isMediumMuon(const reco::Muon&, bool run2016_hip_mitigation=false);
+   bool isSoftMuon(const reco::Muon&, const reco::Vertex&, bool run2016_hip_mitigation=false);
    bool isHighPtMuon(const reco::Muon&, const reco::Vertex&);
+   bool isTrackerHighPtMuon(const reco::Muon&, const reco::Vertex&);
    
    // determine if station was crossed well withing active volume
    unsigned int RequiredStationMask( const reco::Muon& muon,
@@ -116,5 +117,8 @@ namespace muon {
    int sharedSegments( const reco::Muon& muon1, const reco::Muon& muon2, 
                        unsigned int segmentArbitrationMask = reco::MuonSegmentMatch::BestInChamberByDR ) ;
 
+   void setCutBasedSelectorFlags(reco::Muon& muon,
+				 const reco::Vertex* vertex=nullptr,
+				 bool run2016_hip_mitigation=false);
 }
 #endif
diff --git a/DataFormats/MuonReco/src/Muon.cc b/DataFormats/MuonReco/src/Muon.cc
index 295a143025bc1..189e95e64f414 100755
--- a/DataFormats/MuonReco/src/Muon.cc
+++ b/DataFormats/MuonReco/src/Muon.cc
@@ -17,6 +17,7 @@ Muon::Muon(  Charge q, const LorentzVector & p4, const Point & vtx ) :
      type_ = 0;
      bestTunePTrackType_ = reco::Muon::None;
      bestTrackType_ = reco::Muon::None;
+     selectors_ = 0;
 }
 
 Muon::Muon() {
@@ -29,12 +30,12 @@ Muon::Muon() {
    type_ = 0;
    bestTrackType_ = reco::Muon::None;
    bestTunePTrackType_ = reco::Muon::None;
-
+   selectors_ = 0;
 }
 
 bool Muon::overlap( const Candidate & c ) const {
   const RecoCandidate * o = dynamic_cast<const RecoCandidate *>( & c );
-  return ( o != 0 && 
+  return ( o != nullptr && 
 	   ( checkOverlap( track(), o->track() ) ||
 	     checkOverlap( standAloneMuon(), o->standAloneMuon() ) ||
 	     checkOverlap( combinedMuon(), o->combinedMuon() ) ||
@@ -386,8 +387,8 @@ const std::vector<const MuonChamberMatch*> Muon::chambers( int station, int muon
 std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> Muon::pair( const std::vector<const MuonChamberMatch*> &chambers,
      ArbitrationType type ) const
 {
-   MuonChamberMatch* m = 0;
-   MuonSegmentMatch* s = 0;
+   MuonChamberMatch* m = nullptr;
+   MuonSegmentMatch* s = nullptr;
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair(m,s);
 
    if(chambers.empty()) return chamberSegmentPair;
@@ -425,7 +426,7 @@ std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> Muon::pair( const std
 float Muon::dX( int station, int muonSubdetId, ArbitrationType type ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    return chamberSegmentPair.first->x-chamberSegmentPair.second->x;
 }
@@ -434,7 +435,7 @@ float Muon::dY( int station, int muonSubdetId, ArbitrationType type ) const
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    return chamberSegmentPair.first->y-chamberSegmentPair.second->y;
 }
@@ -442,7 +443,7 @@ float Muon::dY( int station, int muonSubdetId, ArbitrationType type ) const
 float Muon::dDxDz( int station, int muonSubdetId, ArbitrationType type ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    return chamberSegmentPair.first->dXdZ-chamberSegmentPair.second->dXdZ;
 }
@@ -451,7 +452,7 @@ float Muon::dDyDz( int station, int muonSubdetId, ArbitrationType type ) const
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    return chamberSegmentPair.first->dYdZ-chamberSegmentPair.second->dYdZ;
 }
@@ -459,7 +460,7 @@ float Muon::dDyDz( int station, int muonSubdetId, ArbitrationType type ) const
 float Muon::pullX( int station, int muonSubdetId, ArbitrationType type, bool includeSegmentError ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    if(includeSegmentError)
       return (chamberSegmentPair.first->x-chamberSegmentPair.second->x)/sqrt(pow(chamberSegmentPair.first->xErr,2)+pow(chamberSegmentPair.second->xErr,2));
@@ -470,7 +471,7 @@ float Muon::pullY( int station, int muonSubdetId, ArbitrationType type, bool inc
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    if(includeSegmentError)
       return (chamberSegmentPair.first->y-chamberSegmentPair.second->y)/sqrt(pow(chamberSegmentPair.first->yErr,2)+pow(chamberSegmentPair.second->yErr,2));
@@ -480,7 +481,7 @@ float Muon::pullY( int station, int muonSubdetId, ArbitrationType type, bool inc
 float Muon::pullDxDz( int station, int muonSubdetId, ArbitrationType type, bool includeSegmentError ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    if(includeSegmentError)
       return (chamberSegmentPair.first->dXdZ-chamberSegmentPair.second->dXdZ)/sqrt(pow(chamberSegmentPair.first->dXdZErr,2)+pow(chamberSegmentPair.second->dXdZErr,2));
@@ -491,7 +492,7 @@ float Muon::pullDyDz( int station, int muonSubdetId, ArbitrationType type, bool
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    if(includeSegmentError)
       return (chamberSegmentPair.first->dYdZ-chamberSegmentPair.second->dYdZ)/sqrt(pow(chamberSegmentPair.first->dYdZErr,2)+pow(chamberSegmentPair.second->dYdZErr,2));
@@ -501,7 +502,7 @@ float Muon::pullDyDz( int station, int muonSubdetId, ArbitrationType type, bool
 float Muon::segmentX( int station, int muonSubdetId, ArbitrationType type ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    return chamberSegmentPair.second->x;
 }
@@ -510,7 +511,7 @@ float Muon::segmentY( int station, int muonSubdetId, ArbitrationType type ) cons
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    return chamberSegmentPair.second->y;
 }
@@ -518,7 +519,7 @@ float Muon::segmentY( int station, int muonSubdetId, ArbitrationType type ) cons
 float Muon::segmentDxDz( int station, int muonSubdetId, ArbitrationType type ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    return chamberSegmentPair.second->dXdZ;
 }
@@ -527,7 +528,7 @@ float Muon::segmentDyDz( int station, int muonSubdetId, ArbitrationType type ) c
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    return chamberSegmentPair.second->dYdZ;
 }
@@ -535,7 +536,7 @@ float Muon::segmentDyDz( int station, int muonSubdetId, ArbitrationType type ) c
 float Muon::segmentXErr( int station, int muonSubdetId, ArbitrationType type ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    return chamberSegmentPair.second->xErr;
 }
@@ -544,7 +545,7 @@ float Muon::segmentYErr( int station, int muonSubdetId, ArbitrationType type ) c
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    return chamberSegmentPair.second->yErr;
 }
@@ -552,7 +553,7 @@ float Muon::segmentYErr( int station, int muonSubdetId, ArbitrationType type ) c
 float Muon::segmentDxDzErr( int station, int muonSubdetId, ArbitrationType type ) const
 {
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasPhi()) return 999999;
    return chamberSegmentPair.second->dXdZErr;
 }
@@ -561,7 +562,7 @@ float Muon::segmentDyDzErr( int station, int muonSubdetId, ArbitrationType type
 {
    if(station==4 && muonSubdetId==MuonSubdetId::DT) return 999999; // no y information
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(chambers(station,muonSubdetId),type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) return 999999;
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) return 999999;
    if(! chamberSegmentPair.second->hasZed()) return 999999;
    return chamberSegmentPair.second->dYdZErr;
 }
@@ -572,7 +573,7 @@ float Muon::trackEdgeX( int station, int muonSubdetId, ArbitrationType type ) co
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -593,7 +594,7 @@ float Muon::trackEdgeY( int station, int muonSubdetId, ArbitrationType type ) co
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -614,7 +615,7 @@ float Muon::trackX( int station, int muonSubdetId, ArbitrationType type ) const
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -635,7 +636,7 @@ float Muon::trackY( int station, int muonSubdetId, ArbitrationType type ) const
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -656,7 +657,7 @@ float Muon::trackDxDz( int station, int muonSubdetId, ArbitrationType type ) con
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -677,7 +678,7 @@ float Muon::trackDyDz( int station, int muonSubdetId, ArbitrationType type ) con
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -698,7 +699,7 @@ float Muon::trackXErr( int station, int muonSubdetId, ArbitrationType type ) con
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -719,7 +720,7 @@ float Muon::trackYErr( int station, int muonSubdetId, ArbitrationType type ) con
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -740,7 +741,7 @@ float Muon::trackDxDzErr( int station, int muonSubdetId, ArbitrationType type )
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -761,7 +762,7 @@ float Muon::trackDyDzErr( int station, int muonSubdetId, ArbitrationType type )
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
@@ -782,7 +783,7 @@ float Muon::trackDist( int station, int muonSubdetId, ArbitrationType type ) con
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
             muonChamber != muonChambers.end(); ++muonChamber) {
@@ -799,7 +800,7 @@ float Muon::trackDistErr( int station, int muonSubdetId, ArbitrationType type )
    if(muonChambers.empty()) return 999999;
 
    std::pair<const MuonChamberMatch*,const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers,type);
-   if(chamberSegmentPair.first==0 || chamberSegmentPair.second==0) {
+   if(chamberSegmentPair.first==nullptr || chamberSegmentPair.second==nullptr) {
       float dist  = 999999;
       float supVar = 999999;
       for(std::vector<const MuonChamberMatch*>::const_iterator muonChamber = muonChambers.begin();
diff --git a/DataFormats/MuonReco/src/MuonSelectors.cc b/DataFormats/MuonReco/src/MuonSelectors.cc
index 08d68fe6273a0..279ff6a86821d 100644
--- a/DataFormats/MuonReco/src/MuonSelectors.cc
+++ b/DataFormats/MuonReco/src/MuonSelectors.cc
@@ -38,7 +38,7 @@ SelectionType selectionTypeFromString( const std::string &label )
       { "ME0MuonArbitrated", ME0MuonArbitrated },
       { "AllGEMMuons", AllGEMMuons },
       { "GEMMuonArbitrated", GEMMuonArbitrated },
-      { 0, (SelectionType)-1 }
+      { nullptr, (SelectionType)-1 }
    };
 
    SelectionType value = (SelectionType)-1;
@@ -361,7 +361,7 @@ bool muon::isGoodMuon( const reco::Muon& muon,
             minNumberOfMatches = 1;
       }
 
-      if(numSegs >= minNumberOfMatches) goodMuon = 1;
+      if(numSegs >= minNumberOfMatches) goodMuon = true;
 
       // Require that last required station have segment
       // If there are zero required stations keep track
@@ -396,11 +396,11 @@ bool muon::isGoodMuon( const reco::Muon& muon,
       detector = lastSegBit < 4 ? 1 : 2;
 
       // Check x information
-      if(fabs(muon.pullX(station,detector,arbitrationType,1)) > maxAbsPullX &&
+      if(fabs(muon.pullX(station,detector,arbitrationType,true)) > maxAbsPullX &&
             fabs(muon.dX(station,detector,arbitrationType)) > maxAbsDx)
          return false;
 
-      if(applyAlsoAngularCuts && fabs(muon.pullDxDz(station,detector,arbitrationType,1)) > maxAbsPullX)
+      if(applyAlsoAngularCuts && fabs(muon.pullDxDz(station,detector,arbitrationType,true)) > maxAbsPullX)
          return false;
 
       // Is this a tight algorithm, i.e. do we bother to check y information?
@@ -408,11 +408,11 @@ bool muon::isGoodMuon( const reco::Muon& muon,
 
          // Check y information
          if (detector == 2) { // CSC
-            if(fabs(muon.pullY(station,2,arbitrationType,1)) > maxAbsPullY &&
+            if(fabs(muon.pullY(station,2,arbitrationType,true)) > maxAbsPullY &&
                   fabs(muon.dY(station,2,arbitrationType)) > maxAbsDy)
                return false;
 
-            if(applyAlsoAngularCuts && fabs(muon.pullDyDz(station,2,arbitrationType,1)) > maxAbsPullY)
+            if(applyAlsoAngularCuts && fabs(muon.pullDyDz(station,2,arbitrationType,true)) > maxAbsPullY)
                return false;
          } else {
             //
@@ -437,12 +437,12 @@ bool muon::isGoodMuon( const reco::Muon& muon,
                if(muon.dY(stationIdx,1,arbitrationType) > 999998) // no y-information
                   continue;
 
-               if(fabs(muon.pullY(stationIdx,1,arbitrationType,1)) > maxAbsPullY &&
+               if(fabs(muon.pullY(stationIdx,1,arbitrationType,true)) > maxAbsPullY &&
                      fabs(muon.dY(stationIdx,1,arbitrationType)) > maxAbsDy) {
                   return false;
                }
 
-               if(applyAlsoAngularCuts && fabs(muon.pullDyDz(stationIdx,1,arbitrationType,1)) > maxAbsPullY)
+               if(applyAlsoAngularCuts && fabs(muon.pullDyDz(stationIdx,1,arbitrationType,true)) > maxAbsPullY)
                   return false;
 
                // If we get this far then great this is a good muon
@@ -485,9 +485,9 @@ bool muon::isGoodMuon( const reco::Muon& muon,
             station  = stationIdx < 4 ? stationIdx+1 : stationIdx-3;
             detector = stationIdx < 4 ? 1 : 2;
 
-            if((fabs(muon.pullX(station,detector,arbitrationType,1)) > maxAbsPullX &&
+            if((fabs(muon.pullX(station,detector,arbitrationType,true)) > maxAbsPullX &&
                   fabs(muon.dX(station,detector,arbitrationType)) > maxAbsDx) ||
-                  (applyAlsoAngularCuts && fabs(muon.pullDxDz(station,detector,arbitrationType,1)) > maxAbsPullX))
+                  (applyAlsoAngularCuts && fabs(muon.pullDxDz(station,detector,arbitrationType,true)) > maxAbsPullX))
                continue;
             else if (detector == 1)
                existsGoodDTSegX = true;
@@ -495,9 +495,9 @@ bool muon::isGoodMuon( const reco::Muon& muon,
             // Is this a tight algorithm?  If yes, use y information
             if (maxAbsDy < 999999) {
                if (detector == 2) { // CSC
-                  if((fabs(muon.pullY(station,2,arbitrationType,1)) > maxAbsPullY &&
+                  if((fabs(muon.pullY(station,2,arbitrationType,true)) > maxAbsPullY &&
                         fabs(muon.dY(station,2,arbitrationType)) > maxAbsDy) ||
-                        (applyAlsoAngularCuts && fabs(muon.pullDyDz(station,2,arbitrationType,1)) > maxAbsPullY))
+                        (applyAlsoAngularCuts && fabs(muon.pullDyDz(station,2,arbitrationType,true)) > maxAbsPullY))
                      continue;
                } else {
 
@@ -506,9 +506,9 @@ bool muon::isGoodMuon( const reco::Muon& muon,
                   else
                      existsDTSegY = true;
 
-                  if((fabs(muon.pullY(station,1,arbitrationType,1)) > maxAbsPullY &&
+                  if((fabs(muon.pullY(station,1,arbitrationType,true)) > maxAbsPullY &&
                         fabs(muon.dY(station,1,arbitrationType)) > maxAbsDy) ||
-                        (applyAlsoAngularCuts && fabs(muon.pullDyDz(station,1,arbitrationType,1)) > maxAbsPullY)) {
+                        (applyAlsoAngularCuts && fabs(muon.pullDyDz(station,1,arbitrationType,true)) > maxAbsPullY)) {
                      continue;
                   }
                }
@@ -856,8 +856,13 @@ bool muon::isLooseMuon(const reco::Muon& muon){
 }
 
 
-bool muon::isMediumMuon(const reco::Muon& muon){
-  if( !( isLooseMuon(muon) && muon.innerTrack()->validFraction() > 0.8 )) return false; 
+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. && 
@@ -867,8 +872,8 @@ bool muon::isMediumMuon(const reco::Muon& muon){
   return (segmentCompatibility(muon) > (goodGlb ? 0.303 : 0.451)); 
 }
 
-
-bool muon::isSoftMuon(const reco::Muon& muon, const reco::Vertex& vtx){
+bool muon::isSoftMuon(const reco::Muon& muon, const reco::Vertex& vtx,
+		      bool run2016_hip_mitigation){
 
   bool muID = muon::isGoodMuon(muon, TMOneStationTight);
 
@@ -881,7 +886,7 @@ bool muon::isSoftMuon(const reco::Muon& muon, const reco::Vertex& vtx){
   
   bool ip = fabs(muon.innerTrack()->dxy(vtx.position())) < 0.3 && fabs(muon.innerTrack()->dz(vtx.position())) < 20.;
   
-  return layers && ip && ishighq;
+  return layers && ip && (ishighq|run2016_hip_mitigation);
 }
 
 
@@ -901,6 +906,21 @@ bool muon::isHighPtMuon(const reco::Muon& muon, const reco::Vertex& vtx){
 
 }
 
+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()/muon.tunePMuonBestTrack()->pt() < 0.3;
+
+  bool ip = fabs(muon.innerTrack()->dxy(vtx.position())) < 0.2 && fabs(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;
    
@@ -929,3 +949,48 @@ int muon::sharedSegments( const reco::Muon& mu, const reco::Muon& mu2, unsigned
   
     return ret; 
 }
+
+void muon::setCutBasedSelectorFlags(reco::Muon& muon, 
+				    const reco::Vertex* 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 dbCorectedRelIso = dbCorrectedIsolation/muon.pt();
+  double tkRelIso = muon.isolationR03().sumPt/muon.pt();
+
+  // Base selectors
+  if (muon::isLooseMuon(muon))        selectors |= reco::Muon::CutBasedIdLoose;
+  if (vertex){
+    if (muon::isTightMuon(muon,*vertex))  selectors |= reco::Muon::CutBasedIdTight;
+    if (muon::isSoftMuon(muon,*vertex,run2016_hip_mitigation))   selectors |= reco::Muon::SoftCutBasedId;
+    if (muon::isHighPtMuon(muon,*vertex)) selectors |= reco::Muon::CutBasedIdGlobalHighPt;
+    if (muon::isTrackerHighPtMuon(muon,*vertex)) selectors |= reco::Muon::CutBasedIdTrkHighPt;
+  }
+  if (muon::isMediumMuon(muon,run2016_hip_mitigation)){
+    selectors |= reco::Muon::CutBasedIdMedium;
+    if ( vertex and 
+	 fabs(muon.muonBestTrack()->dz( vertex->position()))<0.1 and 
+	 fabs(muon.muonBestTrack()->dxy(vertex->position()))< 0.02 )
+      selectors |= reco::Muon::CutBasedIdMediumPrompt;
+  }
+
+  // PF isolation
+  if (dbCorectedRelIso<0.40)    selectors |= reco::Muon::PFIsoVeryLoose;
+  if (dbCorectedRelIso<0.25)    selectors |= reco::Muon::PFIsoLoose;
+  if (dbCorectedRelIso<0.20)    selectors |= reco::Muon::PFIsoMedium;
+  if (dbCorectedRelIso<0.15)    selectors |= reco::Muon::PFIsoTight;
+  if (dbCorectedRelIso<0.10)    selectors |= reco::Muon::PFIsoVeryTight;
+  
+  // Tracker isolation
+  if (tkRelIso<0.10)            selectors |= reco::Muon::TkIsoLoose;
+  if (tkRelIso<0.05)            selectors |= reco::Muon::TkIsoTight;
+
+  muon.setSelectors(selectors);
+}
diff --git a/DataFormats/MuonReco/src/classes_def.xml b/DataFormats/MuonReco/src/classes_def.xml
index eba97ad6aae3f..4d293cdf96636 100644
--- a/DataFormats/MuonReco/src/classes_def.xml
+++ b/DataFormats/MuonReco/src/classes_def.xml
@@ -1,5 +1,6 @@
 <lcgdict>
-  <class name="reco::Muon" ClassVersion="18">
+  <class name="reco::Muon" ClassVersion="19">
+   <version ClassVersion="19" checksum="3201675353"/>
    <version ClassVersion="18" checksum="2257088106"/>
    <version ClassVersion="17" checksum="2296096670"/>
    <version ClassVersion="16" checksum="1808419014"/>
diff --git a/DataFormats/PatCandidates/interface/Muon.h b/DataFormats/PatCandidates/interface/Muon.h
index d72d3d11799f3..803ca2f92668c 100644
--- a/DataFormats/PatCandidates/interface/Muon.h
+++ b/DataFormats/PatCandidates/interface/Muon.h
@@ -265,6 +265,15 @@ namespace pat {
       float pfEcalEnergy() const { return pfEcalEnergy_; }
       void setPfEcalEnergy(float pfEcalEnergy) { pfEcalEnergy_ = pfEcalEnergy; }
 
+      /// near-by jet information
+      float jetPtRatio() const { return jetPtRatio_; }
+      float jetPtRel()   const { return jetPtRel_; }
+      void  setJetPtRatio(float jetPtRatio){ jetPtRatio_ = jetPtRatio; }
+      void  setJetPtRel(float jetPtRel){ jetPtRel_ = jetPtRel; }
+
+      /// Muon MVA
+      float mvaValue() const { return mvaValue_; }
+      void  setMvaValue(float mva){ mvaValue_ = mva; }
     protected:
 
       // ---- for content embedding ----
@@ -332,6 +341,13 @@ namespace pat {
       float puppiNoLeptonsPhotonIso_;
 
       float pfEcalEnergy_;
+
+      /// near-by jet information
+      float jetPtRatio_;
+      float jetPtRel_;
+
+      /// Muon MVA
+      float mvaValue_;
   };
 
 
diff --git a/DataFormats/PatCandidates/src/Muon.cc b/DataFormats/PatCandidates/src/Muon.cc
index 4ea93e6422ce5..f617fc5fd43f4 100644
--- a/DataFormats/PatCandidates/src/Muon.cc
+++ b/DataFormats/PatCandidates/src/Muon.cc
@@ -29,7 +29,10 @@ Muon::Muon() :
     normChi2_(0.0),
     cachedNumberOfValidHits_(false),
     numberOfValidHits_(0),
-    pfEcalEnergy_(0)
+    pfEcalEnergy_(0),
+    jetPtRatio_(0),
+    jetPtRel_(0),
+    mvaValue_(0)
 {
   initImpactParameters();
 }
@@ -53,7 +56,10 @@ Muon::Muon(const reco::Muon & aMuon) :
     normChi2_(0.0),
     cachedNumberOfValidHits_(false),
     numberOfValidHits_(0),
-    pfEcalEnergy_(0)
+    pfEcalEnergy_(0),
+    jetPtRatio_(0),
+    jetPtRel_(0),
+    mvaValue_(0)
 {
   initImpactParameters();
 }
@@ -75,9 +81,12 @@ Muon::Muon(const edm::RefToBase<reco::Muon> & aMuonRef) :
     pfCandidateRef_(),
     cachedNormChi2_(false),
     normChi2_(0.0),
-    cachedNumberOfValidHits_(0),
+    cachedNumberOfValidHits_(false),
     numberOfValidHits_(0),
-    pfEcalEnergy_(0)
+    pfEcalEnergy_(0),
+    jetPtRatio_(0),
+    jetPtRel_(0),
+    mvaValue_(0)
 {
   initImpactParameters();
 }
@@ -99,9 +108,12 @@ Muon::Muon(const edm::Ptr<reco::Muon> & aMuonRef) :
     pfCandidateRef_(),
     cachedNormChi2_(false),
     normChi2_(0.0),
-    cachedNumberOfValidHits_(0),
+    cachedNumberOfValidHits_(false),
     numberOfValidHits_(0),
-    pfEcalEnergy_(0)
+    pfEcalEnergy_(0),
+    jetPtRatio_(0),
+    jetPtRel_(0),
+    mvaValue_(0)
 {
   initImpactParameters();
 }
diff --git a/DataFormats/PatCandidates/src/classes_def_objects.xml b/DataFormats/PatCandidates/src/classes_def_objects.xml
index b4a67c4b4c759..396576d6a9244 100644
--- a/DataFormats/PatCandidates/src/classes_def_objects.xml
+++ b/DataFormats/PatCandidates/src/classes_def_objects.xml
@@ -60,7 +60,8 @@
   </ioread>
 
 
-  <class name="pat::Muon"  ClassVersion="22">
+  <class name="pat::Muon"  ClassVersion="23">
+   <version ClassVersion="23" checksum="3471390119"/>
    <version ClassVersion="22" checksum="1816533594"/>
    <version ClassVersion="21" checksum="1539691612"/>
    <version ClassVersion="20" checksum="357097717"/>
diff --git a/PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h b/PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h
new file mode 100644
index 0000000000000..9a89ab7edd01c
--- /dev/null
+++ b/PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h
@@ -0,0 +1,46 @@
+#ifndef __PhysicsTools_PatAlgos_MuonMvaEstimator__
+#define __PhysicsTools_PatAlgos_MuonMvaEstimator__
+#include "DataFormats/MuonReco/interface/Muon.h"
+#include "DataFormats/PatCandidates/interface/Muon.h"
+#include "DataFormats/BTauReco/interface/JetTag.h"
+#include "TMVA/Reader.h"
+
+namespace reco {
+  class JetCorrector;
+}
+namespace pat {
+  class MuonMvaEstimator{
+  public:
+    MuonMvaEstimator();
+    void initialize(std::string weightsfile, 
+		    float dRmax);
+    void computeMva(const pat::Muon& imuon,
+		    const reco::Vertex& vertex,
+		    const reco::JetTagCollection& bTags,
+		    const reco::JetCorrector* correctorL1=nullptr,
+		    const reco::JetCorrector* correctorL1L2L3Res=nullptr);
+    float mva() const {return mva_;}
+    float jetPtRatio() const {return jetPtRatio_;}
+    float jetPtRel() const {return jetPtRel_;}
+  private:
+    TMVA::Reader tmvaReader_;
+    bool initialized_;
+    float mva_;
+    float dRmax_;
+    
+    /// MVA VAriables
+    float pt_;
+    float eta_;
+    float jetNDauCharged_;
+    float miniRelIsoCharged_;
+    float miniRelIsoNeutral_;
+    float jetPtRel_;
+    float jetPtRatio_;
+    float jetBTagCSV_;
+    float sip_;
+    float log_abs_dxyBS_; 
+    float log_abs_dzPV_;            
+    float segmentCompatibility_;
+  };
+}
+#endif
diff --git a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc
index 74ac80b1855df..1979456865882 100755
--- a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc
+++ b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc
@@ -38,6 +38,9 @@
 #include "FWCore/Utilities/interface/transform.h"
 
 #include "PhysicsTools/PatUtils/interface/MiniIsolation.h"
+#include "PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h"
+#include "FWCore/ParameterSet/interface/FileInPath.h"
+#include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
 
 #include <vector>
 #include <memory>
@@ -47,7 +50,13 @@ using namespace pat;
 using namespace std;
 
 
-PATMuonProducer::PATMuonProducer(const edm::ParameterSet & iConfig) : useUserData_(iConfig.exists("userData")),
+PATMuonProducer::PATMuonProducer(const edm::ParameterSet & iConfig) : 
+  relMiniIsoPUCorrected_(0),
+  useUserData_(iConfig.exists("userData")),
+  computeMuonMVA_(false),
+  recomputeBasicSelectors_(false),
+  mvaDrMax_(0),
+  mvaUseJec_(false),
   isolator_(iConfig.exists("userIsolation") ? iConfig.getParameter<edm::ParameterSet>("userIsolation") : edm::ParameterSet(), consumesCollector(), false)
 {
   // input source
@@ -122,6 +131,7 @@ PATMuonProducer::PATMuonProducer(const edm::ParameterSet & iConfig) : useUserDat
 
   //for mini-isolation calculation
   computeMiniIso_ = iConfig.getParameter<bool>("computeMiniIso");
+
   miniIsoParams_ = iConfig.getParameter<std::vector<double> >("miniIsoParams");
   if(computeMiniIso_ && miniIsoParams_.size() != 9){
       throw cms::Exception("ParameterError") << "miniIsoParams must have exactly 9 elements.\n";
@@ -129,6 +139,27 @@ PATMuonProducer::PATMuonProducer(const edm::ParameterSet & iConfig) : useUserDat
   if(computeMiniIso_)
       pcToken_ = consumes<pat::PackedCandidateCollection >(iConfig.getParameter<edm::InputTag>("pfCandsForMiniIso"));
 
+  // standard selectors
+  recomputeBasicSelectors_ = iConfig.getParameter<bool>("recomputeBasicSelectors");
+  computeMuonMVA_ = iConfig.getParameter<bool>("computeMuonMVA");
+  mvaTrainingFile_ = iConfig.getParameter<std::string>("mvaTrainingFile");
+  if (computeMuonMVA_ and not computeMiniIso_) 
+    throw cms::Exception("ConfigurationError") << "MiniIso is needed for Muon MVA calculation.\n";
+
+  if (computeMuonMVA_) {
+    // pfCombinedInclusiveSecondaryVertexV2BJetTags
+    mvaBTagCollectionTag_  = consumes<reco::JetTagCollection>(iConfig.getParameter<edm::InputTag>("mvaJetTag"));
+    mvaL1Corrector_        = consumes<reco::JetCorrector>(iConfig.getParameter<edm::InputTag>("mvaL1Corrector"));
+    mvaL1L2L3ResCorrector_ = consumes<reco::JetCorrector>(iConfig.getParameter<edm::InputTag>("mvaL1L2L3ResCorrector"));
+    rho_                   = consumes<double>(iConfig.getParameter<edm::InputTag>("rho"));
+    mvaUseJec_ = iConfig.getParameter<bool>("mvaUseJec");
+    mvaDrMax_  = iConfig.getParameter<double>("mvaDrMax");
+
+    // xml training file
+    edm::FileInPath fip(mvaTrainingFile_);
+    mvaEstimator_.initialize(fip.fullPath(),mvaDrMax_);
+  }
+
   // produces vector of muons
   produces<std::vector<Muon> >();
 }
@@ -190,6 +221,16 @@ void PATMuonProducer::produce(edm::Event & iEvent, const edm::EventSetup & iSetu
     iEvent.getByToken(PUPPINoLeptonsIsolation_neutral_hadrons_, PUPPINoLeptonsIsolation_neutral_hadrons);
     iEvent.getByToken(PUPPINoLeptonsIsolation_photons_, PUPPINoLeptonsIsolation_photons);  
   }
+
+  // inputs for muon mva
+  edm::Handle<reco::JetTagCollection> mvaBTagCollectionTag;
+  edm::Handle<reco::JetCorrector> mvaL1Corrector;
+  edm::Handle<reco::JetCorrector> mvaL1L2L3ResCorrector;
+  if (computeMuonMVA_) {
+    iEvent.getByToken(mvaBTagCollectionTag_,mvaBTagCollectionTag);
+    iEvent.getByToken(mvaL1Corrector_,mvaL1Corrector);
+    iEvent.getByToken(mvaL1L2L3ResCorrector_,mvaL1L2L3ResCorrector);
+  }
   
   // prepare the MC genMatchTokens_
   GenAssociations  genMatches(genMatchTokens_.size());
@@ -435,7 +476,64 @@ void PATMuonProducer::produce(edm::Event & iEvent, const edm::EventSetup & iSetu
   // sort muons in pt
   std::sort(patMuons->begin(), patMuons->end(), pTComparator_);
 
-  // put genEvt object in Event
+  // Store standard muon selection decisions and jet related 
+  // quantaties.
+  // Need a separate loop over muons to have all inputs properly
+  // computed and stored in the object.
+  edm::Handle<double> rho;
+  if (computeMuonMVA_) iEvent.getByToken(rho_,rho);
+  const reco::Vertex* pv(nullptr);
+  if (primaryVertexIsValid) pv = &primaryVertex;
+  for(auto& muon: *patMuons){
+    if (recomputeBasicSelectors_){
+      muon.setSelectors(0);
+      bool isRun2016BCDEF = (272728 <= iEvent.run() && iEvent.run() <= 278808);
+      muon::setCutBasedSelectorFlags(muon, pv, isRun2016BCDEF);
+    }
+    if (computeMiniIso_){
+      // MiniIsolation working points
+      double iso = getRelMiniIsoPUCorrected(muon,*rho);
+      muon.setSelector(reco::Muon::MiniIsoLoose,     iso<0.40);
+      muon.setSelector(reco::Muon::MiniIsoMedium,    iso<0.20);
+      muon.setSelector(reco::Muon::MiniIsoTight,     iso<0.10);
+      muon.setSelector(reco::Muon::MiniIsoVeryTight, iso<0.05);
+    }
+    if (computeMuonMVA_ && primaryVertexIsValid){
+      if (mvaUseJec_)
+	mvaEstimator_.computeMva(muon,
+				 primaryVertex,
+				 *(mvaBTagCollectionTag.product()),
+				 &*mvaL1Corrector,
+				 &*mvaL1L2L3ResCorrector);
+      else
+	mvaEstimator_.computeMva(muon,
+				 primaryVertex,
+				 *(mvaBTagCollectionTag.product()));
+      
+      muon.setMvaValue(mvaEstimator_.mva());
+      muon.setJetPtRatio(mvaEstimator_.jetPtRatio());
+      muon.setJetPtRel(mvaEstimator_.jetPtRel());
+      
+      // MVA working points
+      // https://twiki.cern.ch/twiki/bin/viewauth/CMS/LeptonMVA
+      double dB2D  = fabs(muon.dB(pat::Muon::BS2D));
+      double dB3D  = fabs(muon.dB(pat::Muon::PV3D));
+      double edB3D = fabs(muon.edB(pat::Muon::PV3D));
+      double sip3D = edB3D>0?dB3D/edB3D:0.0; 
+      double dz = fabs(muon.muonBestTrack()->dz(primaryVertex.position()));
+
+      // muon preselection
+      if (muon.pt()>5 and muon.isLooseMuon() and
+	  muon.passed(reco::Muon::MiniIsoLoose) and 
+	  sip3D<8.0 and dB2D<0.05 and dz<0.1){
+	muon.setSelector(reco::Muon::MvaLoose,  muon.mvaValue()>-0.60);
+	muon.setSelector(reco::Muon::MvaMedium, muon.mvaValue()>-0.20);
+	muon.setSelector(reco::Muon::MvaTight,  muon.mvaValue()> 0.15);
+      }
+    }
+  }
+
+  // put products in Event
   std::unique_ptr<std::vector<Muon> > ptr(patMuons);
   iEvent.put(std::move(ptr));
 
@@ -528,6 +626,15 @@ void PATMuonProducer::setMuonMiniIso(Muon& aMuon, const PackedCandidateCollectio
   aMuon.setMiniPFIsolation(miniiso);
 }
 
+double PATMuonProducer::getRelMiniIsoPUCorrected(const pat::Muon& muon, float rho)
+{
+  float mindr(miniIsoParams_[0]);
+  float maxdr(miniIsoParams_[1]);
+  float kt_scale(miniIsoParams_[2]);
+  float drcut = pat::miniIsoDr(muon.p4(),mindr,maxdr,kt_scale);
+  return pat::muonRelMiniIsoPUCorrected(muon.miniPFIsolation(), muon.p4(), drcut, rho);
+}
+
 // ParameterSet description for module
 void PATMuonProducer::fillDescriptions(edm::ConfigurationDescriptions & descriptions)
 {
diff --git a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.h b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.h
index ef1c2458c5f66..7b12a93f153f3 100644
--- a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.h
+++ b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.h
@@ -32,7 +32,7 @@
 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
 #include "PhysicsTools/PatAlgos/interface/PATUserDataHelper.h"
 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
-
+#include "PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h"
 
 namespace pat {
   /// foward declarations
@@ -46,9 +46,9 @@ namespace pat {
     /// default constructir
     explicit PATMuonProducer(const edm::ParameterSet & iConfig);
     /// default destructur
-    ~PATMuonProducer();
+    ~PATMuonProducer() override;
     /// everything that needs to be done during the event loop
-    virtual void produce(edm::Event & iEvent, const edm::EventSetup& iSetup) override;
+    void produce(edm::Event & iEvent, const edm::EventSetup& iSetup) override;
     /// description of config file parameters
     static void fillDescriptions(edm::ConfigurationDescriptions & descriptions);
 
@@ -68,7 +68,8 @@ namespace pat {
     /// for the embedding of isoDeposits or userIsolation values
     template<typename T> void readIsolationLabels( const edm::ParameterSet & iConfig, const char* psetName, IsolationLabels& labels, std::vector<edm::EDGetTokenT<edm::ValueMap<T> > > & tokens);
 
-      void setMuonMiniIso(pat::Muon& aMuon, const pat::PackedCandidateCollection *pc);
+    void setMuonMiniIso(pat::Muon& aMuon, const pat::PackedCandidateCollection *pc);
+    double getRelMiniIsoPUCorrected(const pat::Muon& muon, float rho);
 
     // embed various impact parameters with errors
     // embed high level selection
@@ -79,7 +80,8 @@ namespace pat {
 			 bool primaryVertexIsValid,
 			 reco::BeamSpot & beamspot,
 			 bool beamspotIsValid );
-
+    double relMiniIsoPUCorrected( const pat::Muon& aMuon,
+				  double rho);
 
   private:
     /// input source
@@ -89,6 +91,7 @@ namespace pat {
     edm::EDGetTokenT<pat::PackedCandidateCollection > pcToken_;
     bool computeMiniIso_;
     std::vector<double> miniIsoParams_;
+    double relMiniIsoPUCorrected_;
 
     /// embed the track from best muon measurement (global pflow)
     bool embedBestTrack_;
@@ -160,6 +163,17 @@ namespace pat {
     edm::EDGetTokenT<edm::ValueMap<float> > PUPPINoLeptonsIsolation_charged_hadrons_;
     edm::EDGetTokenT<edm::ValueMap<float> > PUPPINoLeptonsIsolation_neutral_hadrons_;
     edm::EDGetTokenT<edm::ValueMap<float> > PUPPINoLeptonsIsolation_photons_;
+    /// standard muon selectors
+    bool computeMuonMVA_;
+    bool recomputeBasicSelectors_;
+    double mvaDrMax_;
+    bool mvaUseJec_;
+    edm::EDGetTokenT<reco::JetTagCollection> mvaBTagCollectionTag_;
+    edm::EDGetTokenT<reco::JetCorrector> mvaL1Corrector_;
+    edm::EDGetTokenT<reco::JetCorrector> mvaL1L2L3ResCorrector_;
+    edm::EDGetTokenT<double> rho_;
+    pat::MuonMvaEstimator mvaEstimator_;
+    std::string mvaTrainingFile_;
     
     /// --- tools ---
     /// comparator for pt ordering
diff --git a/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py b/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py
index be04f46a4458e..b586e677aef3a 100644
--- a/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py
+++ b/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py
@@ -102,6 +102,19 @@
     pfCandsForMiniIso = cms.InputTag("packedPFCandidates"),
     miniIsoParams = cms.vdouble(0.05, 0.2, 10.0, 0.5, 0.0001, 0.01, 0.01, 0.01, 0.0),
 
+    # Standard Muon Selectors and Jet-related observables
+    # Depends on MiniIsolation, so only works in miniaod
+    # Don't forget to set flags properly in miniAOD_tools.py                      
+    computeMuonMVA = cms.bool(False),
+    mvaTrainingFile = cms.string("RecoMuon/MuonIdentification/data/mu_BDTG.weights.xml"),
+    recomputeBasicSelectors = cms.bool(True),
+    mvaUseJec = cms.bool(True),
+    mvaDrMax = cms.double(0.4),
+    mvaJetTag = cms.InputTag("pfCombinedInclusiveSecondaryVertexV2BJetTags"),
+    mvaL1Corrector = cms.InputTag("ak4PFCHSL1FastjetCorrector"),
+    mvaL1L2L3ResCorrector = cms.InputTag("ak4PFCHSL1FastL2L3Corrector"),
+    rho = cms.InputTag("fixedGridRhoFastjetCentralNeutral")
+                          
 )
 
 
diff --git a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py
index 206f6e33f7089..3cb36901b5646 100644
--- a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py
+++ b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py
@@ -26,6 +26,7 @@ def miniAOD_customizeCommon(process):
     process.patMuons.puppiNoLeptonsIsolationPhotons        = cms.InputTag("muonPUPPINoLeptonsIsolation","gamma-DR040-ThresholdVeto000-ConeVeto001")
 
     process.patMuons.computeMiniIso = cms.bool(True)
+    process.patMuons.computeMuonMVA = cms.bool(True)
     
     #
     # disable embedding of electron and photon associated objects already stored by the ReducedEGProducer
diff --git a/PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc b/PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc
new file mode 100644
index 0000000000000..2f827b4f66898
--- /dev/null
+++ b/PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc
@@ -0,0 +1,135 @@
+#include "PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h"
+#include "DataFormats/MuonReco/interface/MuonSelectors.h"
+#include "DataFormats/VertexReco/interface/Vertex.h"
+#include "DataFormats/JetReco/interface/PFJet.h"
+#include "DataFormats/JetReco/interface/PFJetCollection.h"
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
+
+using namespace pat;
+
+MuonMvaEstimator::MuonMvaEstimator():
+  tmvaReader_("!Color:!Silent:Error"),
+  initialized_(false),
+  mva_(0),
+  dRmax_(0)
+{}
+
+void MuonMvaEstimator::initialize(std::string weightsfile,
+				  float dRmax)
+{ 
+  if (initialized_) return;
+  tmvaReader_.AddVariable("LepGood_pt",                    &pt_               );
+  tmvaReader_.AddVariable("LepGood_eta",                   &eta_              );
+  tmvaReader_.AddVariable("LepGood_jetNDauChargedMVASel",  &jetNDauCharged_   );
+  tmvaReader_.AddVariable("LepGood_miniRelIsoCharged",     &miniRelIsoCharged_);
+  tmvaReader_.AddVariable("LepGood_miniRelIsoNeutral",     &miniRelIsoNeutral_);
+  tmvaReader_.AddVariable("LepGood_jetPtRelv2",            &jetPtRel_         );
+  tmvaReader_.AddVariable("min(LepGood_jetPtRatiov2,1.5)", &jetPtRatio_       );
+  tmvaReader_.AddVariable("max(LepGood_jetBTagCSV,0)",     &jetBTagCSV_       );
+  tmvaReader_.AddVariable("LepGood_sip3d",                 &sip_              );
+  tmvaReader_.AddVariable("log(abs(LepGood_dxy))",         &log_abs_dxyBS_    ); 
+  tmvaReader_.AddVariable("log(abs(LepGood_dz))",          &log_abs_dzPV_     );
+  tmvaReader_.AddVariable("LepGood_segmentCompatibility",  &segmentCompatibility_);
+  tmvaReader_.BookMVA("BDTG",weightsfile);
+  dRmax_ = dRmax;
+  initialized_ = true;
+};
+
+float ptRel(const reco::Candidate::LorentzVector& muP4, 
+	    const reco::Candidate::LorentzVector& jetP4, 
+	    bool subtractMuon=true) 
+{
+  reco::Candidate::LorentzVector jp4 = jetP4;
+  if (subtractMuon) jp4-=muP4;
+  float dot = muP4.Vect().Dot( jp4.Vect() );
+  float ptrel = muP4.P2() - dot*dot/jp4.P2();
+  ptrel = ptrel>0 ? sqrt(ptrel) : 0.0;
+  return ptrel;
+}
+
+void MuonMvaEstimator::computeMva(const pat::Muon& muon,
+				  const reco::Vertex& vertex,
+				  const reco::JetTagCollection& bTags,
+				  const reco::JetCorrector* correctorL1,
+				  const reco::JetCorrector* correctorL1L2L3Res)
+{
+  if (not initialized_) 
+    throw cms::Exception("FatalError") << "MuonMVA is not initialized";
+  pt_                   = muon.pt();
+  eta_                  = muon.eta();
+  segmentCompatibility_ = muon.segmentCompatibility();
+  miniRelIsoCharged_ = muon.miniPFIsolation().chargedHadronIso();
+  miniRelIsoNeutral_ = muon.miniPFIsolation().neutralHadronIso();
+
+  double dB2D  = fabs(muon.dB(pat::Muon::BS2D));
+  double dB3D  = muon.dB(pat::Muon::PV3D);
+  double edB3D = muon.edB(pat::Muon::PV3D);
+  double dz    = fabs(muon.muonBestTrack()->dz(vertex.position()));
+  sip_  = edB3D>0?fabs(dB3D/edB3D):0.0; 
+  log_abs_dxyBS_     = dB2D>0?log(dB2D):0; 
+  log_abs_dzPV_      = dz>0?log(dz):0;
+
+  //Initialise loop variables
+  double minDr = 9999;
+  double jecL1L2L3Res = 1.;
+  double jecL1 = 1.;
+
+  jetPtRatio_ = -99;
+  jetPtRel_   = -99;
+  jetBTagCSV_ = -999;
+  jetNDauCharged_ = -1;
+
+  for (const auto& tagI: bTags){
+    // for each muon with the lepton 
+    double dr = deltaR(*(tagI.first), muon);
+    if(dr > minDr) continue;  
+    minDr = dr;
+      
+    const reco::Candidate::LorentzVector& muP4(muon.p4()); 
+    reco::Candidate::LorentzVector jetP4(tagI.first->p4());
+
+    if (correctorL1 && correctorL1L2L3Res){
+      jecL1L2L3Res = correctorL1L2L3Res->correction(*(tagI.first));
+      jecL1 = correctorL1->correction(*(tagI.first));
+    }
+
+    // Get b-jet info
+    jetBTagCSV_ = tagI.second;
+    jetNDauCharged_ = 0;
+    for (auto jet: tagI.first->getJetConstituentsQuick()){
+      const reco::PFCandidate *pfcand = dynamic_cast<const reco::PFCandidate*>(jet);
+      if (pfcand==nullptr) throw cms::Exception("ConfigurationError") << "Cannot get jet constituents";
+      if (pfcand->charge()==0) continue;
+      auto bestTrackPtr = pfcand->bestTrack();
+      if (!bestTrackPtr) continue;
+      if (!bestTrackPtr->quality(reco::Track::highPurity)) continue;
+      if (bestTrackPtr->pt()<1.) continue;
+      if (bestTrackPtr->hitPattern().numberOfValidHits()<8) continue;
+      if (bestTrackPtr->hitPattern().numberOfValidPixelHits()<2) continue;
+      if (bestTrackPtr->normalizedChi2()>=5) continue;
+
+      if (std::fabs(bestTrackPtr->dxy(vertex.position())) > 0.2) continue;
+      if (std::fabs(bestTrackPtr->dz(vertex.position())) > 17) continue;
+      jetNDauCharged_++;
+    }
+
+    if(minDr < dRmax_){
+      if ((jetP4-muP4).Rho()<0.0001){ 
+	jetPtRel_ = 0;
+	jetPtRatio_ = 1;
+      } else {
+	jetP4 -= muP4/jecL1;
+	jetP4 *= jecL1L2L3Res;
+	jetP4 += muP4;
+      
+	jetPtRatio_ = muP4.pt()/jetP4.pt();
+	jetPtRel_ = ptRel(muP4,jetP4);
+      }
+    }
+  }
+
+  if (jetPtRatio_>1.5) jetPtRatio_ = 1.5;
+  if (jetBTagCSV_<0) jetBTagCSV_ = 0;
+  mva_ = tmvaReader_.EvaluateMVA("BDTG");
+};
diff --git a/PhysicsTools/PatUtils/interface/MiniIsolation.h b/PhysicsTools/PatUtils/interface/MiniIsolation.h
index 3bfd1080d47b8..3d14ab8275f60 100644
--- a/PhysicsTools/PatUtils/interface/MiniIsolation.h
+++ b/PhysicsTools/PatUtils/interface/MiniIsolation.h
@@ -15,13 +15,20 @@
 
 namespace pat{
 
-    // see src file for definitions of parameters
-    PFIsolation getMiniPFIsolation(const pat::PackedCandidateCollection *pfcands, const math::XYZTLorentzVector& p4,
+  float miniIsoDr(const math::XYZTLorentzVector &p4, float mindr, float maxdr,
+		  float kt_scale);
+
+  // see src file for definitions of parameters
+  PFIsolation getMiniPFIsolation(const pat::PackedCandidateCollection *pfcands, const math::XYZTLorentzVector& p4,
                                    float mindr=0.05, float maxdr=0.2, float kt_scale=10.0,
                                    float ptthresh=0.5, float deadcone_ch=0.0001,
                                    float deadcone_pu=0.01, float deadcone_ph=0.01, float deadcone_nh=0.01,
                                    float dZ_cut=0.0);
 
+  float muonRelMiniIsoPUCorrected(const PFIsolation& iso,
+				  const math::XYZTLorentzVector& p4,
+				  float dr,
+				  float rho);
 }
 
 #endif
diff --git a/PhysicsTools/PatUtils/src/MiniIsolation.cc b/PhysicsTools/PatUtils/src/MiniIsolation.cc
index d47313d3998d6..a0f575528e2cd 100644
--- a/PhysicsTools/PatUtils/src/MiniIsolation.cc
+++ b/PhysicsTools/PatUtils/src/MiniIsolation.cc
@@ -11,6 +11,12 @@ namespace pat{
 // Excludes PFCands inside of "deadcone" radius.
 // For nh, ph, pu, only include particles with pT > ptthresh
 // Some documentation can be found here: https://hypernews.cern.ch/HyperNews/CMS/get/susy/1991.html
+
+float miniIsoDr(const math::XYZTLorentzVector &p4, float mindr, float maxdr,
+		 float kt_scale){
+  return std::max(mindr, std::min(maxdr, float(kt_scale/p4.pt())));
+}
+
 PFIsolation getMiniPFIsolation(const pat::PackedCandidateCollection *pfcands,
                                  const math::XYZTLorentzVector &p4, float mindr, float maxdr,
                                  float kt_scale, float ptthresh, float deadcone_ch,
@@ -19,7 +25,7 @@ PFIsolation getMiniPFIsolation(const pat::PackedCandidateCollection *pfcands,
 {
     
     float chiso=0, nhiso=0, phiso=0, puiso=0;
-    float drcut = std::max(mindr, std::min(maxdr, float(kt_scale/p4.pt())));
+    float drcut = miniIsoDr(p4,mindr,maxdr,kt_scale);
     for(auto const & pc : *pfcands){
         float dr = deltaR(p4, pc.p4());
         if(dr>drcut)
@@ -49,4 +55,22 @@ PFIsolation getMiniPFIsolation(const pat::PackedCandidateCollection *pfcands,
 
 }
 
+  float muonRelMiniIsoPUCorrected(const PFIsolation& iso,
+				  const math::XYZTLorentzVector& p4,
+				  float dr,
+				  float rho)
+  {
+    double absEta = fabs(p4.eta());
+    double ea = 0;
+    //Spring15 version
+    if      (absEta<0.800) ea = 0.0735;
+    else if (absEta<1.300) ea = 0.0619;
+    else if (absEta<2.000) ea = 0.0465;
+    else if (absEta<2.200) ea = 0.0433;
+    else if (absEta<2.500) ea = 0.0577;
+    double correction = rho * ea * (dr/0.3) * (dr/0.3);
+    double correctedIso = iso.chargedHadronIso() + std::max(0.0, iso.neutralHadronIso()+iso.photonIso() - correction);
+    return correctedIso/p4.pt();
+  }
+
 }
diff --git a/RecoHI/HiMuonAlgos/python/HiRecoMuon_cff.py b/RecoHI/HiMuonAlgos/python/HiRecoMuon_cff.py
index 6d70f7ea0e5bc..a0e331db7ca70 100644
--- a/RecoHI/HiMuonAlgos/python/HiRecoMuon_cff.py
+++ b/RecoHI/HiMuonAlgos/python/HiRecoMuon_cff.py
@@ -46,6 +46,7 @@
 muons.FillSelectorMaps = cms.bool(False)
 muons.FillShoweringInfo = cms.bool(False)
 muons.FillCosmicsIdMap = cms.bool(False)
+muons.vertices = cms.InputTag("hiSelectedVertex")
 muonRecoHighLevelPbPb = cms.Sequence(muons)
 
 # HI muon sequence (passed to RecoHI.Configuration.Reconstruction_HI_cff)
diff --git a/RecoMuon/MuonIdentification/plugins/MuonProducer.cc b/RecoMuon/MuonIdentification/plugins/MuonProducer.cc
index c93903adeb5fb..076ad9d33f79d 100644
--- a/RecoMuon/MuonIdentification/plugins/MuonProducer.cc
+++ b/RecoMuon/MuonIdentification/plugins/MuonProducer.cc
@@ -7,6 +7,9 @@
 #include "RecoMuon/MuonIdentification/plugins/MuonProducer.h"
 
 #include "RecoMuon/MuonIsolation/interface/MuPFIsoHelper.h"
+#include "DataFormats/VertexReco/interface/VertexFwd.h"
+#include "DataFormats/VertexReco/interface/Vertex.h"
+#include "DataFormats/MuonReco/interface/MuonSelectors.h"
 
 
 #include <boost/foreach.hpp>
@@ -48,6 +51,7 @@ MuonProducer::MuonProducer(const edm::ParameterSet& pSet):debug_(pSet.getUntrack
   fillDetectorBasedIsolation_ = pSet.getParameter<bool>("FillDetectorBasedIsolation"); 
   fillShoweringInfo_          = pSet.getParameter<bool>("FillShoweringInfo");
   fillTimingInfo_             = pSet.getParameter<bool>("FillTimingInfo");
+  computeStandardSelectors_   = pSet.getParameter<bool>("ComputeStandardSelectors");
 
   produces<reco::MuonCollection>();
 
@@ -172,6 +176,11 @@ MuonProducer::MuonProducer(const edm::ParameterSet& pSet):debug_(pSet.getUntrack
     }
     
   }
+  
+  if (computeStandardSelectors_){
+    vertexes_ = consumes<reco::VertexCollection>(pSet.getParameter<edm::InputTag>("vertices"));
+  }
+
 }
 
 /// Destructor
@@ -196,6 +205,13 @@ void MuonProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup)
    edm::Handle<reco::PFCandidateCollection> pfCandidates; 
    event.getByToken(thePFCandToken_, pfCandidates);
 
+   edm::Handle<reco::VertexCollection> primaryVertices;
+   const reco::Vertex* vertex(nullptr);
+   if (computeStandardSelectors_){
+     event.getByToken(vertexes_, primaryVertices);
+     if (!primaryVertices->empty())
+       vertex = &(primaryVertices->front());
+   }
 
    // fetch collections for PFIso
    if(fillPFIsolation_) thePFIsoHelper->beginEvent(event);
@@ -426,6 +442,13 @@ void MuonProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup)
        cosmicIdColl[i] = (*cosmicIdMap)[muRef];
        cosmicCompColl[i] = (*cosmicCompMap)[muRef];
      }
+     
+     // Standard Selectors - keep it at the end so that all inputs are available
+     if (computeStandardSelectors_){
+       outMuon.setSelectors(0); // reset flags
+       bool isRun2016BCDEF = (272728 <= event.run() && event.run() <= 278808);
+       muon::setCutBasedSelectorFlags(outMuon, vertex, isRun2016BCDEF);
+     }
 
      outputMuons->push_back(outMuon); 
      ++i;
diff --git a/RecoMuon/MuonIdentification/plugins/MuonProducer.h b/RecoMuon/MuonIdentification/plugins/MuonProducer.h
index 4db86985e79e3..314df96d9be86 100644
--- a/RecoMuon/MuonIdentification/plugins/MuonProducer.h
+++ b/RecoMuon/MuonIdentification/plugins/MuonProducer.h
@@ -39,6 +39,7 @@ namespace reco {class Track;}
 #include "DataFormats/Common/interface/ValueMap.h"
 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
+#include "DataFormats/VertexReco/interface/VertexFwd.h"
 
 
 class MuPFIsoHelper;
@@ -51,10 +52,10 @@ class MuonProducer : public edm::stream::EDProducer<> {
   MuonProducer(const edm::ParameterSet&);
 
   /// Destructor
-  virtual ~MuonProducer();
+  ~MuonProducer() override;
 
   /// reconstruct muons
-  virtual void produce(edm::Event&, const edm::EventSetup&) override;
+  void produce(edm::Event&, const edm::EventSetup&) override;
 
 
   typedef std::vector<edm::InputTag> InputTags;
@@ -98,6 +99,7 @@ class MuonProducer : public edm::stream::EDProducer<> {
   bool fillDetectorBasedIsolation_;
   bool fillShoweringInfo_;
   bool fillTimingInfo_;
+  bool computeStandardSelectors_;
 
   edm::InputTag theTrackDepositName;
   edm::InputTag theEcalDepositName;
@@ -135,6 +137,8 @@ class MuonProducer : public edm::stream::EDProducer<> {
 
   std::vector<std::map<std::string,edm::InputTag> > pfIsoMapNames;
   std::vector<std::map<std::string,edm::EDGetTokenT<edm::ValueMap<double> > > > pfIsoMapTokens_;
+
+  edm::EDGetTokenT<reco::VertexCollection> vertexes_;
   
 };
 #endif
diff --git a/RecoMuon/MuonIdentification/python/muons_cfi.py b/RecoMuon/MuonIdentification/python/muons_cfi.py
index 2d4fc80f7aaa5..2decc0fed6117 100644
--- a/RecoMuon/MuonIdentification/python/muons_cfi.py
+++ b/RecoMuon/MuonIdentification/python/muons_cfi.py
@@ -91,5 +91,9 @@
                        ShowerInfoMap = cms.InputTag("muonShowerInformation"),
 
                        FillCosmicsIdMap = cms.bool(True),
-                       CosmicIdMap = cms.InputTag("cosmicsVeto")
+                       CosmicIdMap = cms.InputTag("cosmicsVeto"),
+
+                       ComputeStandardSelectors = cms.bool(True),
+                       vertices = cms.InputTag("offlinePrimaryVertices")
+                       
                        )