diff --git a/PhysicsAnalysis/MuonID/MuonSelectorTools/MuonSelectorTools/MuonSelectionTool.h b/PhysicsAnalysis/MuonID/MuonSelectorTools/MuonSelectorTools/MuonSelectionTool.h
index fcbf2551f3fd9222291403f9dbcd4c9c751b138b..cfa746580d54e6d0b8c5ba2e9f38b5cf41035039 100644
--- a/PhysicsAnalysis/MuonID/MuonSelectorTools/MuonSelectorTools/MuonSelectionTool.h
+++ b/PhysicsAnalysis/MuonID/MuonSelectorTools/MuonSelectorTools/MuonSelectionTool.h
@@ -119,8 +119,8 @@ namespace CP {
         /// This is necessary only when the resolution is very optimistic in the MC such that a large smearing is applied
         bool passedBMVmimicCut(const xAOD::Muon&) const;
 
-	///Returns a sorted vector of chamber indices
-	std::vector<int> getChamberIndicesSorted(const xAOD::Muon& mu) const;
+	///Returns a vector of the muon's segments, sorted according to chamber index
+	std::vector<const xAOD::MuonSegment*> getSegmentsSorted(const xAOD::Muon& mu) const;
 
         /// Maximum pseudorapidity for the selected muons
         double m_maxEta;
@@ -139,6 +139,7 @@ namespace CP {
         bool m_useAllAuthors;
         bool m_use2stationMuonsHighPt;
         bool m_useMVALowPt;
+	bool m_useSegmentTaggedLowPt;
         bool m_doBadMuonVetoMimic;
 
         SG::ReadHandleKey<xAOD::EventInfo> m_eventInfo{this, "EventInfoContName", "EventInfo", "event info key"};
@@ -148,6 +149,10 @@ namespace CP {
         std::string m_MVAreaderFile_EVEN_MuGirl;
         std::string m_MVAreaderFile_ODD_MuGirl;
 
+	std::string m_MVAreaderFile_MuTagIMO_etaBin1;
+	std::string m_MVAreaderFile_MuTagIMO_etaBin2;
+	std::string m_MVAreaderFile_MuTagIMO_etaBin3;
+
         std::string m_BMVcutFile;
 
         /// Checks for each histogram
@@ -182,9 +187,10 @@ namespace CP {
         std::unique_ptr<TMVA::Reader> m_readerO_MUID;
         std::unique_ptr<TMVA::Reader> m_readerE_MUGIRL;
         std::unique_ptr<TMVA::Reader> m_readerO_MUGIRL;
-
-        // TMVA initialize function
-        void PrepareReader(TMVA::Reader* reader);
+	
+	std::unique_ptr<TMVA::Reader> m_reader_MUTAGIMO_etaBin1;
+	std::unique_ptr<TMVA::Reader> m_reader_MUTAGIMO_etaBin2;
+	std::unique_ptr<TMVA::Reader> m_reader_MUTAGIMO_etaBin3;
 
         // variables for the TMVA readers
         mutable std::mutex m_low_pt_mva_mutex;
diff --git a/PhysicsAnalysis/MuonID/MuonSelectorTools/Root/MuonSelectionTool.cxx b/PhysicsAnalysis/MuonID/MuonSelectorTools/Root/MuonSelectionTool.cxx
index 60100ea19c0eeb2abd607e82c350b2382127e3ac..b8f018d19b658cb9d73b79da655b708026985ec1 100644
--- a/PhysicsAnalysis/MuonID/MuonSelectorTools/Root/MuonSelectionTool.cxx
+++ b/PhysicsAnalysis/MuonID/MuonSelectorTools/Root/MuonSelectionTool.cxx
@@ -39,12 +39,12 @@ namespace {
       return chamberIndexOrder;
     }
 
-    //This is the comparison function for the sorting of chamber indices
-    bool chamberIndexCompare(int first, int second) {
+    //This is the comparison function for the sorting of segments according to the chamber index
+  bool chamberIndexCompare(const xAOD::MuonSegment* first, const xAOD::MuonSegment* second) {
 
       static const std::vector<int> chamberIndexOrder = initializeChamberIdxOrder();
       
-      return (chamberIndexOrder[first] < chamberIndexOrder[second]);
+      return (chamberIndexOrder[first->chamberIndex()] < chamberIndexOrder[second->chamberIndex()]);
     }
 
     static const SG::AuxElement::Accessor<float> mePt_acc("MuonSpectrometerPt");
@@ -71,8 +71,9 @@ namespace CP {
         // for users of high-pT working point to choose whether to include "safe" 2-station muons
         declareProperty("Use2stationMuonsHighPt", m_use2stationMuonsHighPt = true);
 
-        // for users of low-pT working point to choose whether to use MVA
+        // for users of low-pT working point to choose whether to use MVA and whether to include MuTagIMO muons
         declareProperty("UseMVALowPt", m_useMVALowPt = false);
+	declareProperty( "UseSegmentTaggedLowPt", m_useSegmentTaggedLowPt = false );
 
         // MVA configs for low-pT working point
         declareProperty(
@@ -87,6 +88,12 @@ namespace CP {
         declareProperty(
             "MVAreaderFile_ODD_MuGirl",
             m_MVAreaderFile_ODD_MuGirl = "MuonSelectorTools/190118_PrelimLowPtMVA/LowPtMVA_Weights/BDTG_9JAN2019_MuGirl_ODD.weights.xml");
+	declareProperty( "MVAreaderFile_MuTagIMO_etaBin1", m_MVAreaderFile_MuTagIMO_etaBin1 = 
+			 "dev/MuonSelectorTools/181121_MuTagIMO_BDT/BDT_NOV2021_MuTagIMO_etaBin1.weights.xml");
+	declareProperty( "MVAreaderFile_MuTagIMO_etaBin2", m_MVAreaderFile_MuTagIMO_etaBin2 = 
+			 "dev/MuonSelectorTools/181121_MuTagIMO_BDT/BDT_NOV2021_MuTagIMO_etaBin2.weights.xml");
+	declareProperty( "MVAreaderFile_MuTagIMO_etaBin3", m_MVAreaderFile_MuTagIMO_etaBin3 = 
+			 "dev/MuonSelectorTools/181121_MuTagIMO_BDT/BDT_NOV2021_MuTagIMO_etaBin3.weights.xml");
 
         // switch to cut away the tail of very large smearing in MC to mimic the effect of the bad muon veto for 2-station muons in the
         // high-pT selection
@@ -227,6 +234,33 @@ namespace CP {
             m_readerE_MUGIRL = make_mva_reader(weightPath_EVEN_MuGirl);
 
             m_readerO_MUGIRL = make_mva_reader(weightPath_ODD_MuGirl);
+
+	    if (m_useSegmentTaggedLowPt) {
+
+	      TString weightPath_MuTagIMO_etaBin1 = PathResolverFindCalibFile(m_MVAreaderFile_MuTagIMO_etaBin1);
+	      TString weightPath_MuTagIMO_etaBin2 = PathResolverFindCalibFile(m_MVAreaderFile_MuTagIMO_etaBin2);
+	      TString weightPath_MuTagIMO_etaBin3 = PathResolverFindCalibFile(m_MVAreaderFile_MuTagIMO_etaBin3);
+
+	      auto make_mva_reader_MuTagIMO = [](TString file_path, bool useSeg2ChamberIndex) {
+
+                std::vector<std::string> mva_var_names;
+		if (useSeg2ChamberIndex) mva_var_names.push_back("muonSeg2ChamberIndex");
+		mva_var_names.push_back("muonSeg1ChamberIndex");
+		mva_var_names.push_back("muonSeg1NPrecisionHits");
+		mva_var_names.push_back("muonSegmentDeltaEta");
+		mva_var_names.push_back("muonSeg1GlobalR");
+		mva_var_names.push_back("muonSeg1Chi2OverDoF");
+		mva_var_names.push_back("muonSCS");
+		
+                std::unique_ptr<TMVA::Reader> reader = std::make_unique<TMVA::Reader>(mva_var_names);
+                reader->BookMVA("BDT", file_path);
+                return reader;
+	      };
+
+	      m_reader_MUTAGIMO_etaBin1 = make_mva_reader_MuTagIMO(weightPath_MuTagIMO_etaBin1, false);
+	      m_reader_MUTAGIMO_etaBin2 = make_mva_reader_MuTagIMO(weightPath_MuTagIMO_etaBin2, false);
+	      m_reader_MUTAGIMO_etaBin3 = make_mva_reader_MuTagIMO(weightPath_MuTagIMO_etaBin3, true);
+	    }
         }
         ATH_CHECK(m_eventInfo.initialize());
 
@@ -663,15 +697,33 @@ namespace CP {
             return false;
         }
 
-        // requiring combined muons
-        if (mu.muonType() != xAOD::Muon::Combined) {
+        // requiring combined muons, unless segment-tags are included
+	if (!m_useSegmentTaggedLowPt) {
+	  if (mu.muonType() != xAOD::Muon::Combined) {
             ATH_MSG_VERBOSE("Muon is not combined - fail low-pT");
             return false;
-        }
-        if (mu.author() != xAOD::Muon::MuGirl && mu.author() != xAOD::Muon::MuidCo) {
+	  }
+	}
+	else {
+	  if (mu.muonType() != xAOD::Muon::Combined && mu.muonType() != xAOD::Muon::SegmentTagged) {
+            ATH_MSG_VERBOSE("Muon is not combined or segment-tagged - fail low-pT");
+            return false;
+	  }
+	}
+
+	// author check
+	if (!m_useSegmentTaggedLowPt) {
+	  if (mu.author() != xAOD::Muon::MuGirl && mu.author() != xAOD::Muon::MuidCo) {
             ATH_MSG_VERBOSE("Muon is neither MuGirl nor MuidCo - fail low-pT");
             return false;
-        }
+	  }
+	}
+	else {
+	  if (mu.author() != xAOD::Muon::MuGirl && mu.author() != xAOD::Muon::MuidCo && mu.author() != xAOD::Muon::MuTagIMO) {
+            ATH_MSG_VERBOSE("Muon is neither MuGirl / MuidCo / MuTagIMO - fail low-pT");
+            return false;
+	  }
+	}
 
         // applying Medium selection above pT = 18 GeV
         if (mu.pt() * MeVtoGeV > 18.) {
@@ -703,17 +755,19 @@ namespace CP {
         }
 
         // requiring explicitely >=1 station (2 in the |eta|>1.3 region when Medium selection is not explicitely required)
-        uint8_t nprecisionLayers{0};
-        if (!mu.summaryValue(nprecisionLayers, xAOD::SummaryType::numberOfPrecisionLayers)) {
+	if( mu.muonType() == xAOD::Muon::Combined) {
+	  uint8_t nprecisionLayers{0};
+	  if (!mu.summaryValue(nprecisionLayers, xAOD::SummaryType::numberOfPrecisionLayers)) {
             ATH_MSG_WARNING("passedLowPtEfficiencyCuts - #precision layers missing in combined muon! Failing selection...");
             return false;
-        }
-        uint nStationsCut = (std::abs(mu.eta()) > 1.3 && std::abs(mu.eta()) < 1.55) ? 2 : 1;
-        if (nprecisionLayers < nStationsCut) {
+	  }
+	  uint nStationsCut = (std::abs(mu.eta()) > 1.3 && std::abs(mu.eta()) < 1.55) ? 2 : 1;
+	  if (nprecisionLayers < nStationsCut) {
             ATH_MSG_VERBOSE("number of precision layers = " << (int)nprecisionLayers << " is lower than cut value " << nStationsCut
-                                                            << " - fail low-pT");
+			    << " - fail low-pT");
             return false;
-        }
+	  }
+	}
 
         // reject MuGirl muon if not found also by MuTagIMO
         if (m_useAllAuthors) {
@@ -758,22 +812,22 @@ namespace CP {
     }
 
 
-    std::vector<int> MuonSelectionTool::getChamberIndicesSorted(const xAOD::Muon& mu) const {
+  std::vector<const xAOD::MuonSegment*> MuonSelectionTool::getSegmentsSorted(const xAOD::Muon& mu) const {
 
-      std::vector<int> chamberIndices_sorted;
-      chamberIndices_sorted.reserve(mu.nMuonSegments());
+      std::vector<const xAOD::MuonSegment*> segments_sorted;
+      segments_sorted.reserve(mu.nMuonSegments());
       
       for (unsigned int i = 0; i < mu.nMuonSegments(); i++) {
 	
 	if (mu.muonSegment(i) == nullptr)
 	  ATH_MSG_WARNING("The muon reports more segments than are available. Please report this to the muon software community!");
 	else
-	  chamberIndices_sorted.push_back(mu.muonSegment(i)->chamberIndex());
+	  segments_sorted.push_back(mu.muonSegment(i));
       }
 
-      sort(chamberIndices_sorted.begin(), chamberIndices_sorted.end(), chamberIndexCompare);
+      sort(segments_sorted.begin(), segments_sorted.end(), chamberIndexCompare);
 
-      return chamberIndices_sorted;
+      return segments_sorted;
     }
 
 
@@ -796,27 +850,66 @@ namespace CP {
             return false;
         }
 
-        float segChamberIdx1{-1.}, segChamberIdx2{-1}, middleHoles{-1};
+        float seg1ChamberIdx{-1.}, seg2ChamberIdx{-1.}, middleHoles{-1.}, seg1NPrecisionHits{-1.}, seg1GlobalR{-1.}, seg1Chi2OverDoF{-1.};
+
+	std::vector<const xAOD::MuonSegment*> muonSegments = getSegmentsSorted(mu);
+
+	if (mu.author() == xAOD::Muon::MuTagIMO && muonSegments.size() == 0)
+	  ATH_MSG_WARNING("passedLowPtEfficiencyMVACut - found segment-tagged muon with no segments!");
 
-	std::vector<int> chamberIndices = getChamberIndicesSorted(mu);
+        seg1ChamberIdx =     (muonSegments.size() > 0) ? muonSegments[0]->chamberIndex() : -9;
+        seg2ChamberIdx =     (muonSegments.size() > 1) ? muonSegments[1]->chamberIndex() : -9;
+
+	//these variables are only used for MuTagIMO
+	if (mu.author() == xAOD::Muon::MuTagIMO) {
+	  seg1NPrecisionHits = (muonSegments.size() > 0) ? muonSegments[0]->nPrecisionHits() : -1;
+	  seg1GlobalR =        (muonSegments.size() > 0) ? sqrt( muonSegments[0]->x()*muonSegments[0]->x() + 
+								 muonSegments[0]->y()*muonSegments[0]->y() +
+								 muonSegments[0]->z()*muonSegments[0]->z() ) : 0;
+	  seg1Chi2OverDoF =    (muonSegments.size() > 0) ? muonSegments[0]->chiSquared() / muonSegments[0]->numberDoF() : -1;
+	}
 
-        segChamberIdx1 = (chamberIndices.size() > 0) ? chamberIndices[0] : -9;
-        segChamberIdx2 = (chamberIndices.size() > 1) ? chamberIndices[1] : -9;
         middleHoles = middleSmallHoles + middleLargeHoles;
 
         // get event number from event info
         SG::ReadHandle<xAOD::EventInfo> eventInfo(m_eventInfo);
 
-        const std::vector<float> var_vector{
-            momentumBalanceSig, 
-            CurvatureSig, 
-            scatteringNeigbour, 
-            energyLoss, 
-            middleHoles, 
-            muonSegmentDeltaEta, 
-            segChamberIdx1, 
-            segChamberIdx2,
-        };
+	// variables for the BDT
+        std::vector<float> var_vector;
+	if (mu.author() == xAOD::Muon::MuidCo || mu.author() == xAOD::Muon::MuGirl) {
+	  var_vector = {
+	    momentumBalanceSig, 
+	    CurvatureSig, 
+	    scatteringNeigbour, 
+	    energyLoss, 
+	    middleHoles, 
+	    muonSegmentDeltaEta, 
+	    seg1ChamberIdx, 
+	    seg2ChamberIdx
+	  };
+	}
+	else {
+	  if ( std::abs(mu.eta()) >= 1.3 )
+	    var_vector = {
+	      seg2ChamberIdx,
+	      seg1ChamberIdx,
+	      seg1NPrecisionHits,
+	      muonSegmentDeltaEta,
+	      seg1GlobalR,
+	      seg1Chi2OverDoF,
+	      std::abs(CurvatureSig)
+	    };
+	  else
+	    var_vector = {
+	      seg1ChamberIdx,
+	      seg1NPrecisionHits,
+	      muonSegmentDeltaEta,
+	      seg1GlobalR,
+	      seg1Chi2OverDoF,
+	      std::abs(CurvatureSig)
+	    };
+	}
+
         // use different trainings for even/odd numbered events
         TMVA::Reader *reader_MUID, *reader_MUGIRL;
         if (eventInfo->eventNumber() % 2 == 1) {
@@ -827,20 +920,33 @@ namespace CP {
             reader_MUGIRL = m_readerO_MUGIRL.get();
         }
 
+	//BDT for MuTagIMO is binned in |eta|
+	TMVA::Reader *reader_MUTAGIMO;
+	if ( std::abs(mu.eta()) < 0.7 )
+	  reader_MUTAGIMO = m_reader_MUTAGIMO_etaBin1.get();
+	else if ( std::abs(mu.eta()) < 1.3 )
+	  reader_MUTAGIMO = m_reader_MUTAGIMO_etaBin2.get();
+	else
+	  reader_MUTAGIMO = m_reader_MUTAGIMO_etaBin3.get();
+
         // get the BDT discriminant response
         float BDTdiscriminant;
 
         if (mu.author() == xAOD::Muon::MuidCo)
             BDTdiscriminant = reader_MUID->EvaluateMVA(var_vector, "BDTG");
-        else if (mu.author() == xAOD::Muon::MuGirl || mu.author() == xAOD::Muon::MuTagIMO)
+        else if (mu.author() == xAOD::Muon::MuGirl)
             BDTdiscriminant = reader_MUGIRL->EvaluateMVA(var_vector, "BDTG");
+	else if (mu.author() == xAOD::Muon::MuTagIMO && m_useSegmentTaggedLowPt)
+	    BDTdiscriminant = reader_MUTAGIMO->EvaluateMVA(var_vector, "BDT");
         else {
             ATH_MSG_WARNING("Invalid author for low-pT MVA, failing selection...");
             return false;
         }
 
-        // cut on dicriminant at -0.6
-        if (BDTdiscriminant > -0.6) {
+	//cut on dicriminant
+	float BDTcut = (mu.author() == xAOD::Muon::MuTagIMO) ? 0.12 : -0.6;
+
+        if (BDTdiscriminant > BDTcut) {
             ATH_MSG_VERBOSE("Passed low-pT MVA cut");
             return true;
         } else {
@@ -1176,7 +1282,7 @@ namespace CP {
         if (mu.muonType() == xAOD::Muon::CaloTagged && std::abs(mu.eta()) < 0.105)
             return true;  // removed the passedCaloTagQuality(mu) until this is better understood in r22
         // ::
-        if (mu.muonType() == xAOD::Muon::SegmentTagged && std::abs(mu.eta()) < 0.105) return true;
+        if (mu.muonType() == xAOD::Muon::SegmentTagged && (std::abs(mu.eta()) < 0.105 || m_useSegmentTaggedLowPt)) return true;
         // ::
         if (mu.author() == xAOD::Muon::MuidSA && std::abs(mu.eta()) > 2.4) return true;
         // ::