diff --git a/Muon/MuonTrackRec/src/component/MuonCombRec.cpp b/Muon/MuonTrackRec/src/component/MuonCombRec.cpp
index 59d5821f3f7c6ef06ca01385e5b48ed714f5538b..2ff18db7bc6b4ed94329791779298ef476ec5f9b 100644
--- a/Muon/MuonTrackRec/src/component/MuonCombRec.cpp
+++ b/Muon/MuonTrackRec/src/component/MuonCombRec.cpp
@@ -38,7 +38,7 @@ DECLARE_COMPONENT( MuonCombRec )
 
 namespace {
   // -- initialize the pad size. Hardwired to speed up.
-  constexpr std::array<double, 20> m_padX = {{  
+  constexpr std::array<double, 20> m_padX = {{
   //   R1         R2         R3         R4
     10.1667 ,  20.1667 ,  40.333  ,  80.6667 ,   // M1
      6.41667,  12.75   ,  25.5    ,  51.     ,   // M2
@@ -47,7 +47,7 @@ namespace {
     31.6667 ,  63.     , 126.     , 252      }}; // M5
 
   //             R1         R2         R3         R4
-  constexpr std::array<double, 20> m_padY = {{ 
+  constexpr std::array<double, 20> m_padY = {{
   //   R1         R2         R3         R4
     24.7586 ,  49.7584 ,  99.7581 , 199.757  ,   // M1
     31.3585 ,  62.9583 , 126.158  , 252.557  ,   // M2
@@ -157,9 +157,18 @@ StatusCode MuonCombRec::initialize() {
   // incident service
   incSvc()->addListener( this, IncidentType::EndEvent );
 
+  // OL 171114: offset introduced
+  // the idea is that the various arrays will always be indexed following the M1, M2 etc. enums
+
   // geometry stuff --->
   m_nStation = m_muonDetector->stations();
-  if(m_nStation != 5){
+  if ( m_nStation == 5 ) {
+    // Run 1 or 2
+    m_stationOffset = 0;
+  } else if ( m_nStation == 4 ) {
+    // Upgrade
+    m_stationOffset = 1;
+  } else {
     err()<<"Wrong station number: "<<m_nStation<<endmsg;
     return StatusCode::FAILURE;
   }
@@ -173,9 +182,9 @@ StatusCode MuonCombRec::initialize() {
   }
 
   // -- get the z position of stations
-  for (int station=0; station<m_nStation; ++station){
-    m_zStations[station] = m_muonDetector->getStationZ(station);
-    if ( msgLevel(MSG::DEBUG) ) debug() << "Z of station M" << station+1 << ": " << m_zStations[station] << endmsg;
+  for (int station = 0; station < m_nStation; ++station){
+    m_zStations[station + m_stationOffset] = m_muonDetector->getStationZ(station);
+    if ( msgLevel(MSG::DEBUG) ) debug() << "Z of station M" << station+m_stationOffset+1 << ": " << m_zStations[station+m_stationOffset] << endmsg;
   }
   // <---
 
@@ -196,15 +205,18 @@ StatusCode MuonCombRec::initialize() {
   // local stuff
   if ( msgLevel(MSG::DEBUG) ){
     debug()<<"xFOIs: "<<endmsg;
-    for(int is = 0; is < m_nStation-1; is++){
+    for(int is = m_stationOffset; is < m_nStation + m_stationOffset - 1; is++) { // [M1] M2-4
+      debug() << "M" << is+1 << " ";
       for(int ir = 0; ir < m_nRegion; ir++){
         int key = is*4+ir;
         debug()<< m_xFOIs[key] <<" ";
       }
       debug()<<endmsg;
     }
+
     debug()<<"yFOIs: "<<endmsg;
-    for(int is = 0; is < m_nStation-1; is++){
+    for(int is = m_stationOffset; is < m_nStation + m_stationOffset - 1; is++) {
+      debug() << "M" << is+1 << " ";
       for(int ir = 0; ir < m_nRegion; ir++){
         int key = is*4+ir;
         debug()<< m_yFOIs[key] <<" ";
@@ -215,7 +227,8 @@ StatusCode MuonCombRec::initialize() {
 
   if ( msgLevel(MSG::DEBUG) ){
     debug()<< "x pad sizes: " << endmsg;
-    for(int is = 0; is < m_nStation; is++){
+    for(int is = m_stationOffset; is < m_nStation + m_stationOffset; is++){
+      debug() << "M" << is+1 << " ";
       for(int ir = 0; ir < m_nRegion; ir++){
         int key = is*4+ir;
         debug()<< m_padX[key] << " ";
@@ -223,7 +236,8 @@ StatusCode MuonCombRec::initialize() {
       debug() << endmsg;
     }
     debug()<< "y pad sizes: " << endmsg;
-    for(int is = 0; is < m_nStation; is++){
+    for(int is = m_stationOffset; is < m_nStation + m_stationOffset; is++){
+      debug() << "M" << is+1 << " ";
       for(int ir = 0; ir < m_nRegion; ir++){
         int key = is*4+ir;
         debug()<< m_padY[key] << " ";
@@ -389,9 +403,9 @@ StatusCode MuonCombRec::muonSearch() {
   // -- Local array x-talk hits
   std::vector<const MuonHit*> xtHits;
 
-  int SS = m_seedStation; // SS = Seed Station
+  int SS = m_seedStation; // SS = Seed Station (numbering scheme where 0 = M1, 1 = M2 etc.
 
-  if ( msgLevel(MSG::DEBUG) ) debug() << "Seed Station is: " << SS << endmsg;
+  if ( msgLevel(MSG::DEBUG) ) debug() << "Seed Station is: M" << SS+1 << endmsg;
 
   // -- Loop over SS hits and copy them in a local array with SS seeds
   for ( int region = 0; region < m_nRegion; ++region){
@@ -408,7 +422,7 @@ StatusCode MuonCombRec::muonSearch() {
     candidates[SS].clear();
     candidates[SS].push_back(seedHit);
 
-    if ( msgLevel(MSG::DEBUG) ) debug() << "Seed : " << seedHit->station() << " "
+    if ( msgLevel(MSG::DEBUG) ) debug() << "Seed : M" << seedHit->station() + m_stationOffset + 1 << " "
 					<< seedHit->region() << " "
 					<< seedHit->X() << " "
 					<< seedHit->Y() << " "
@@ -429,7 +443,7 @@ StatusCode MuonCombRec::muonSearch() {
       });
       if(used!=xtHits.end()) {
           if ( msgLevel(MSG::DEBUG) ) debug() << " Hit already used for another seed "
-                                              << (*used)->station() << " "
+                                              << "M" << ((*used)->station() + m_stationOffset + 1) << " "
                                               << (*used)->region() << " "
                                               << (*used)->X() << " "
                                               << (*used)->Y() << " "
@@ -451,11 +465,11 @@ StatusCode MuonCombRec::muonSearch() {
 
         const MuonHit* seedHit2 = *itSSn;
 
-        int istam = seedHit->station();
+        int istam = seedHit->station() + m_stationOffset; // in upgrade station() returns 0 for M2
         int iregm = seedHit->region();
         int keym = istam*4+iregm;
 
-        int istan = seedHit2->station();
+        int istan = seedHit2->station() + m_stationOffset;
         int iregn = seedHit2->region();
         int keyn = istan*4+iregn;
 
@@ -481,7 +495,7 @@ StatusCode MuonCombRec::muonSearch() {
         debug() << "Accumulated x-talk points for this seed: -->" << endmsg;
 
         for( const MuonHit* hit : seedXTHits ){
-          debug() << hit->station() << " " << hit->region()
+          debug() << "M" << hit->station() + m_stationOffset + 1 << " " << hit->region()
                   << " " << hit->X()
                   << " " << hit->Y()
                   << " " << hit->Z() <<endmsg;
@@ -549,7 +563,7 @@ StatusCode MuonCombRec::muonSearch() {
 
     bool matchingFound = true;
 
-    for(int is = nextStation; is >= 0; --is){
+    for(int is = nextStation; is >= m_stationOffset; --is){
 
       if(is == m_skipStation) continue; // skip a station
 
@@ -646,7 +660,7 @@ StatusCode MuonCombRec::muonSearch() {
 // match hits with the track extrapolation
 //========================================================================
 StatusCode MuonCombRec::findMatching( const double x, const double y,
-                                      const int station, const int regionBefore,
+                                      const int station, const int regionBefore, // station follows the M1=0, M2=1 scheme
                                       std::vector<const MuonHit*> &candidates ){
 
   double toll = m_tolForRegion[regionBefore];
@@ -716,11 +730,11 @@ StatusCode MuonCombRec::findMatching( const double x, const double y,
 
         if(deltaX<2&&deltaY<2) continue; // same pad...;-)
 
-        int csta = candidate->station();
+        int csta = candidate->station() + m_stationOffset;
         int creg = candidate->region();
         int ckey = csta*4+creg;
 
-        int ista = hit->station();
+        int ista = hit->station() + m_stationOffset;
         int ireg = hit->region();
         int ikey = ista*4+ireg;
 
@@ -748,7 +762,7 @@ StatusCode MuonCombRec::cloneKiller()
 
       for(const auto& hitT1 : hitsT1) { // loop over hits of 1st trk
 
-        int st1 = hitT1->station();
+        int st1 = hitT1->station() + m_stationOffset;
         if(st1 != M2 && st1 != M3 && st1 != M4) continue;//look only in M2,M3,M4
 
         for(const auto& hitT2 : hitsT2){ // loop over hits of 2nd trk
@@ -806,7 +820,7 @@ StatusCode MuonCombRec::strongCloneKiller()
 
       for(const auto& hitT1 : hitsT1){ // loop over hits of 1st trk
 
-        int st1 = hitT1->station();
+        int st1 = hitT1->station() + m_stationOffset;
         if(st1 != M2 && st1 != M3) continue; // look only in M2,M3
 
         for(const auto& hitT2 : hitsT2){ // loop over hits of 2nd trk
@@ -855,7 +869,7 @@ StatusCode MuonCombRec::sortMuonHits(){
 
   // -- Fill the sorted hit array for later track reconstruction
   for( const auto& hit : *m_trackhits ){
-    int station = hit.station();
+    int station = hit.station() + m_stationOffset;
     int region  = hit.region();
     m_sortedHits[station*4+region].push_back(&hit);
   }
@@ -863,7 +877,7 @@ StatusCode MuonCombRec::sortMuonHits(){
   if ( msgLevel(MSG::DEBUG) ){
     for( const auto& hitVec : m_sortedHits){
       for( const MuonHit* hit : hitVec ){
-        debug() << "sorted hit stored station: " << hit->station()
+        debug() << "sorted hit stored station: M" << hit->station() + m_stationOffset + 1
                 <<" sorted hit stored  region: "<< hit->region()
                 <<" x,y,z: "
                 << hit->X() << " "
diff --git a/Muon/MuonTrackRec/src/component/MuonCombRec.h b/Muon/MuonTrackRec/src/component/MuonCombRec.h
index 5600f80aa7051c974354e86eba702f4a8f525959..0d44f766c4d9ff67c306cb5c3d96cf2d50c35a7b 100644
--- a/Muon/MuonTrackRec/src/component/MuonCombRec.h
+++ b/Muon/MuonTrackRec/src/component/MuonCombRec.h
@@ -141,7 +141,8 @@ private:
 
   void clear();
 
-  int m_nStation;
+  int m_stationOffset; // offset into lookup tables (0 in Run1/2, 1 in upgrade)
+  int m_nStation; // number of stations (5 in Run1/2, 4 in upgrade)
   int m_nRegion;
 
   // -- timing
@@ -156,48 +157,48 @@ private:
 
   // enable timers
   Gaudi::Property<bool> m_measureTime {this, "MeasureTime", false};
-  
+
   // enable clone finding and killing
   bool m_cloneKiller;
-  Gaudi::Property<bool> m_optCloneKiller 
-    {this, "CloneKiller", true, 
+  Gaudi::Property<bool> m_optCloneKiller
+    {this, "CloneKiller", true,
     "enable clone finding and killing"};
 
   // enable strong clone finding and killing
   bool m_strongCloneKiller;
-  Gaudi::Property<bool> m_optStrongCloneKiller 
+  Gaudi::Property<bool> m_optStrongCloneKiller
     {this, "StrongCloneKiller", true,
     "enable strong clone finding and killing"};
 
   // station to be skipped (e.g. for eff. study): 0-4;
   // -1 = keep all stns (default)
   int m_skipStation;
-  Gaudi::Property<int> m_optSkipStation 
+  Gaudi::Property<int> m_optSkipStation
     {this, "SkipStation", -1,
     "station to be skipped (e.g. for eff. study): 0-4; -1 = keep all stns (default)"};
 
-  Gaudi::Property<bool> m_physicsTiming 
-    {this, "PhysicsTiming", true, 
+  Gaudi::Property<bool> m_physicsTiming
+    {this, "PhysicsTiming", true,
     "if true, we assume that timing is the final one (accounting for TOF from primary vx)"};
 
-  Gaudi::Property<bool> m_assumeCosmics 
-    {this, "AssumeCosmics", false, 
+  Gaudi::Property<bool> m_assumeCosmics
+    {this, "AssumeCosmics", false,
     "if true (default) we assume that tracks are of cosmic origin (can be backward)"};
 
-  Gaudi::Property<bool> m_assumePhysics 
-    {this, "AssumePhysics", true, 
+  Gaudi::Property<bool> m_assumePhysics
+    {this, "AssumePhysics", true,
     "if true we assume that tracks have the 'right' direction (pz>0)"};
 
   Gaudi::Property<std::string> m_decToolName
     {this, "DecodingTool", "MuonHitDecode"
     "name of decoding tool (MuonHitDecode for offline, MuonMonHitDecode for online monitoring)"};
 
-  Gaudi::Property<std::string> m_padToolName 
-    {this, "PadRecTool", "MuonPadRec", 
+  Gaudi::Property<std::string> m_padToolName
+    {this, "PadRecTool", "MuonPadRec",
     "name of pad rec tool (MuonPadRec only option so far)"};
-  
-  Gaudi::Property<std::string> m_clusterToolName 
-    {this, "ClusterTool", "MuonFakeClustering", 
+
+  Gaudi::Property<std::string> m_clusterToolName
+    {this, "ClusterTool", "MuonFakeClustering",
     "name of clustering tool"};
 
   /// cross talk
@@ -207,8 +208,8 @@ private:
   int m_seedStation;
   Gaudi::Property<int> m_optSeedStation {this, "SeedStation", 4, "default seed: M5"};
 
-  Gaudi::Property<std::string> m_trackOutputLoc 
-    {this, "TracksOutputLocation", "Rec/Muon/Track", 
+  Gaudi::Property<std::string> m_trackOutputLoc
+    {this, "TracksOutputLocation", "Rec/Muon/Track",
     "LHCb tracks output location in TES"};
 
   /// main steering reconstruction routine
diff --git a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp
index 9ad18771522a7b828b08a55fc92f1bd924c7855a..9b6d93fb6aeb4349830b9de1ac20f7b57d66bae1 100644
--- a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp
+++ b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp
@@ -50,14 +50,34 @@ void MuonPadFromCoord::handle ( const Incident& incident )
   }
 }
 
-
-
-
 //=============================================================================
 StatusCode 	MuonPadFromCoord::initialize ()
 {
   StatusCode sc = GaudiTool::initialize() ;
   if (!sc) return sc;
+
+  // -- retrieve the pointer to the muon detector
+  if ( msgLevel(MSG::DEBUG) ) debug() << "Retrieve MuonDet" << endmsg;
+  m_muonDetector = getDet<DeMuonDetector>(DeMuonLocation::Default);
+  if(!m_muonDetector){
+    error()<<"error retrieving Muon Detector geometry "<<endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  // geometry stuff --->
+  auto nStations = m_muonDetector->stations();
+  if ( nStations == 5 ) {
+    // Run 1 or 2
+    m_stationOffset = 0;
+  } else if ( nStations == 4 ) {
+    // Upgrade
+    m_stationOffset = 1;
+    debug() << "Got " << nStations << " muon stations => setting offset to 1" << endmsg;
+  } else {
+    err() << "Wrong station number: " << nStations << endmsg;
+    return StatusCode::FAILURE;
+  }
+
   incSvc()->addListener( this, IncidentType::EndEvent );
   return StatusCode::SUCCESS ;
 }
@@ -89,23 +109,23 @@ StatusCode MuonPadFromCoord::buildLogicalPads(const std::vector<MuonLogHit>& ) {
   m_pads.reserve( 2*coords->size() );
   // loop over the coords
   for (const auto&  coord : *coords) {
-    MuonTileID tile = MuonTileID(coord->key());
-    unsigned int station = tile.station();
-    unsigned int region = tile.region();
-    double uncross = (station == 0 || ((station>2)&&(region==0))) ? false : coord->uncrossed();
+    auto tile = MuonTileID(coord->key());
+    auto station = tile.station() + m_stationOffset; // offset means 0 means M1 even in the upgrade
+    auto region = tile.region();
+    auto uncross = (station == 0 || ((station>2)&&(region==0))) ? false : coord->uncrossed();
     int time1= (int) coord->digitTDC1();
     int time2= (int) coord->digitTDC2();
 
     if(uncross) {
-      m_tiles.push_back(tile);
+      m_tiles.push_back( tile );
       m_hits.emplace_back( tile );
-      m_hits.back().setTime(time1);
-      m_pads.emplace_back(&m_hits.back());  // m_hits was reserved with enough capacity, so yes, this works
+      m_hits.back().setTime( time1 );
+      m_pads.emplace_back( &m_hits.back() );  // m_hits was reserved with enough capacity, so yes, this works
     } else {
       const std::vector<MuonTileID > thetiles= coord->digitTile();
       tile = MuonTileID(thetiles[0]);
-      m_tiles.push_back(tile);
-      m_hits.emplace_back(tile);
+      m_tiles.push_back( tile );
+      m_hits.emplace_back( tile );
       auto hit1 = &m_hits.back(); // m_hits was reserved with enough capacity, so yes, this works
       m_hits.back().setTime(time1);
       if (thetiles.size() == 1) {
@@ -125,5 +145,3 @@ StatusCode MuonPadFromCoord::buildLogicalPads(const std::vector<MuonLogHit>& ) {
   m_padsReconstructed=true;
   return StatusCode::SUCCESS;
 }
-
-
diff --git a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h
index 88f0923435a36b00151133ad5ad2825cf77b6aa7..a803cd3b96a8ca60927541e6a969d2bb5d7bd887 100644
--- a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h
+++ b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h
@@ -1,4 +1,4 @@
-#ifndef MUONPADFROMCOORD_H 
+#ifndef MUONPADFROMCOORD_H
 #define MUONPADFROMCOORD_H 1
 #include <map>
 #include "GaudiAlg/GaudiTool.h"
@@ -6,7 +6,7 @@
 #include "MuonInterfaces/IMuonPadRec.h"            // Interface
 
 /** @class MuonPadFromCoord MuonPadFromCoord.h
- *  converts MuonCoord to MuonLogPad objects 
+ *  converts MuonCoord to MuonLogPad objects
  *
  *  @author G.Graziani
  *  @date   2014-03-11
@@ -19,9 +19,9 @@ namespace LHCB {
 
 class MuonPadFromCoord final : public extends<GaudiTool, IMuonPadRec , IIncidentListener>
 {
-public: 
+public:
   /// Standard constructor
-  MuonPadFromCoord( const std::string& type, 
+  MuonPadFromCoord( const std::string& type,
               const std::string& name,
               const IInterface* parent);
 
@@ -35,11 +35,13 @@ public:
   // from IIncidentListener
   void handle ( const Incident& incident ) override;
 private:
-  
+
   std::vector<LHCb::MuonTileID> m_tiles;
   std::vector<MuonLogHit> m_hits;
   std::vector<MuonLogPad> m_pads;
   void clearPads();
+  unsigned int m_stationOffset = 0;
   bool m_padsReconstructed = false;
+  DeMuonDetector* m_muonDetector = nullptr;
 };
 #endif // MUONPADFROMCOORD_H
diff --git a/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp b/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp
index 83e78cf0e70202633a310ab4eaaf05532d09cc5f..08cb6e45932fe9af231ae6fb77ef6643195d8f64 100644
--- a/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp
+++ b/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp
@@ -14,6 +14,7 @@
 
 // from TrackEvent
 #include "Event/TrackParameters.h"
+#include "Event/StateParameters.h"
 #include "Event/State.h"
 #include "Event/TrackTypes.h"
 #include "Event/Track.h"
@@ -71,7 +72,7 @@ StatusCode  MuonTrackMomRec::initBdl() {
   // compute integrated B field before M1 at x=y=0
   // we will use this value for all tracks
   Gaudi::XYZPoint  begin( 0., 0., 0. );
-  Gaudi::XYZPoint  end( 0., 0.,  m_muonDetector->getStationZ(0) );
+  Gaudi::XYZPoint  end( 0., 0., StateParameters::ZEndRich2 );
   Gaudi::XYZVector bdl;
 
   m_bIntegrator -> calculateBdlAndCenter(begin, end, 0.0001,
@@ -89,7 +90,7 @@ StatusCode  MuonTrackMomRec::initBdl() {
 std::unique_ptr<LHCb::Track> MuonTrackMomRec::recMomentum(MuonTrack& track) const
 {
   StatusCode sc =  StatusCode::SUCCESS;
-  double Zfirst = m_muonDetector->getStationZ(0);
+  const auto Zfirst = StateParameters::ZEndRich2;
   // create a state at the Z of M1
   Gaudi::XYZPoint trackPos(track.bx() + track.sx() * Zfirst,
                            track.by() + track.sy() * Zfirst,
diff --git a/Pr/PrMCTools/src/PrMuonStubAssociator.cpp b/Pr/PrMCTools/src/PrMuonStubAssociator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..59ba6f8dbe4ce6dd04c943c32c8cdbf8a41a38d6
--- /dev/null
+++ b/Pr/PrMCTools/src/PrMuonStubAssociator.cpp
@@ -0,0 +1,108 @@
+// Include files
+
+
+// from Gaudi
+#include "GaudiKernel/IDataManagerSvc.h"
+#include "GaudiKernel/SmartIF.h"
+
+// from Event
+#include "Event/Track.h"
+#include "Linker/LinkerWithKey.h"
+
+#include "Event/MCParticle.h"
+#include "Event/MCVertex.h"
+
+// local
+#include "PrMuonStubAssociator.h"
+
+DECLARE_ALGORITHM_FACTORY( PrMuonStubAssociator )
+
+//=============================================================================
+// Standard constructor, initializes variables
+//=============================================================================
+  PrMuonStubAssociator::PrMuonStubAssociator( const std::string& name, ISvcLocator* pSvcLocator )
+  : GaudiAlgorithm( name, pSvcLocator )
+{
+}
+
+//=============================================================================
+// Destructor
+//=============================================================================
+PrMuonStubAssociator::~PrMuonStubAssociator() { }
+
+//=============================================================================
+// Initialisation. Check parameters
+//=============================================================================
+StatusCode PrMuonStubAssociator::initialize()
+{
+  // Mandatory initialization of GaudiAlgorithm
+  StatusCode sc = GaudiAlgorithm::initialize();
+  if( sc.isFailure() ) {
+    return sc;
+  }
+
+  info() << "Processing container: " << m_container.value() << endmsg;
+
+  m_muonPad2MC = tool<IMuonPad2MCTool>("MuonPad2MCTool", this);
+  if ( !m_muonPad2MC ) {
+    error() << "Failed to create muonPad2MC tool!" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+//=============================================================================
+// Main execution
+//=============================================================================
+StatusCode PrMuonStubAssociator::execute()
+{
+  // Retrieve the MCParticles
+  auto mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default );
+  if ( !mcParts ) {
+    if( msgLevel(MSG::ERROR) ) {
+      error() << "MCParticles at " << LHCb::MCParticleLocation::Default << "' do not exist" <<endmsg;
+    }
+    return StatusCode::SUCCESS;
+  }
+
+  // Create the Linker table from Track to MCParticle
+  // Sorted by decreasing weight, so first retrieved has highest weight
+  // This has to be done, even if there are no tracks in the event, to satisfy the DST writer
+  LinkerWithKey<LHCb::MCParticle,LHCb::Track> myLinker( evtSvc(), msgSvc(), m_container );
+  myLinker.reset();
+
+  // Retrieve the Tracks
+  auto tracks = getIfExists<LHCb::Track::Range>( m_container );
+  if ( tracks.empty() ) {
+    if( msgLevel(MSG::DEBUG) ) {
+      debug() << "No tracks in container '" << m_container.value() << "'. Skipping." <<endmsg;
+    }
+    return StatusCode::SUCCESS;
+  } else if( msgLevel(MSG::DEBUG) ) {
+    debug() << "Loaded " << tracks.size() << " tracks from '" << m_container.value() << "'" << endmsg;
+  }
+
+  for ( const auto& track : tracks ) {
+    auto nMuon = 0;
+    std::map<LHCb::MCParticle*,unsigned int> matches;
+    for ( const auto& id : track->lhcbIDs() ) {
+      if ( id.isMuon() ) {
+        nMuon++;
+        auto mc_part = m_muonPad2MC->Pad2MC( id.muonID() );
+        if ( mc_part && mc_part->parent() == mcParts ) {
+          matches[ mc_part ]++;
+        }
+      }
+    }
+
+    for ( const auto& match : matches ) {
+      auto fraction = 1.0 * match.second / nMuon;
+      if ( fraction > m_matchThreshold ) {
+        myLinker.link( track, match.first, fraction );
+      }
+    }
+  }
+
+  return StatusCode::SUCCESS;
+}
diff --git a/Pr/PrMCTools/src/PrMuonStubAssociator.h b/Pr/PrMCTools/src/PrMuonStubAssociator.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ac429d6a7387866c5b865b1aeb4b63763b6ae95
--- /dev/null
+++ b/Pr/PrMCTools/src/PrMuonStubAssociator.h
@@ -0,0 +1,41 @@
+#ifndef PRMUONSTUBASSOCIATOR_H
+#define PRMUONSTUBASSOCIATOR_H 1
+
+// Include files
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "MCInterfaces/IMuonPad2MCTool.h"
+
+/** @class TrackAssociator TrackAssociator.h
+ *
+ *  This algorithm computes the link between a Track and a MCParticle.
+ *  The requirement is a match of both the Velo/VP and the T part of the
+ *  Track. If there are not enough coordinates, the match is assumed so that
+ *  a Velo only or a T only are matched properly.
+ *  The required fraction of hits is a jobOption 'FractionOK', default 0.70.
+ *
+ *  Rewritten for the upgrade, handles all containers in one instance
+ *
+ *  @author Olivier Callot
+ *  @date   2012-04-04
+ */
+
+class PrMuonStubAssociator : public GaudiAlgorithm {
+public:
+
+  // Standard constructor
+  PrMuonStubAssociator( const std::string& name, ISvcLocator* pSvcLocator );
+
+  // Destructor
+  virtual ~PrMuonStubAssociator();
+
+  StatusCode initialize() override;  ///< Algorithm initialization
+  StatusCode execute() override;     ///< Algorithm execution
+
+private:
+  Gaudi::Property<std::string> m_container { this, "Container", "Rec/Track/StandaloneMuonFitted" };
+  Gaudi::Property<float> m_matchThreshold { this, "MatchThreshold", 0.7 };
+
+  IMuonPad2MCTool* m_muonPad2MC = nullptr;
+};
+
+#endif // PRMUONSTUBASSOCIATOR_H
diff --git a/Tr/TrackTools/src/MuonSeeding.cpp b/Tr/TrackTools/src/MuonSeeding.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5eac7293983fa75cbca20a1672143293ebebb539
--- /dev/null
+++ b/Tr/TrackTools/src/MuonSeeding.cpp
@@ -0,0 +1,184 @@
+// Include files
+
+// local
+#include "MuonSeeding.h"
+#include "Event/VPCluster.h"
+
+//-----------------------------------------------------------------------------
+// Implementation file for class : MuonSeeding
+//
+// 2010-09-14 : Michel De Cian
+//-----------------------------------------------------------------------------
+
+// Declaration of the Algorithm Factory
+DECLARE_ALGORITHM_FACTORY( MuonSeeding )
+
+MuonSeeding::MuonSeeding( const std::string& name, ISvcLocator* pSvcLocator )
+  : GaudiAlgorithm( name, pSvcLocator )
+{
+  declareProperty( "MuonRecTool", m_trackTool );
+  declareProperty( "Fitter", m_trackFitter );
+  declareProperty( "MuomMomTool", m_momentumTool );
+  declareProperty( "Extrapolator", m_extrapolator );
+}
+
+//=============================================================================
+// Initialization
+//=============================================================================
+StatusCode MuonSeeding::initialize() {
+  const StatusCode sc = GaudiAlgorithm::initialize(); // must be executed first
+  if ( sc.isFailure() ) {
+    return sc;  // error printed already by GaudiAlgorithm
+  }
+
+  if ( msgLevel(MSG::DEBUG) ) {
+    debug() << "==> Initialize" << endmsg;
+  }
+
+  // -- Tracking tools
+  m_trackTool.retrieve();
+  m_extrapolator.retrieve();
+  m_trackFitter.retrieve();
+  m_momentumTool.retrieve();
+
+  return sc;
+}
+
+//=============================================================================
+// Main execution
+//=============================================================================
+StatusCode MuonSeeding::execute()
+{
+  //Chrono chrono( chronoSvc(),name()+"::execute()" );
+
+  if ( msgLevel(MSG::DEBUG) ) debug() << "==> Execute" << endmsg;
+
+  // -- This is where we assume the track came from
+  std::vector<double> PVPos = { 0.0, 0.0, 0.0 }; // TODO take the VELO position?
+
+  /*const LHCb::VPCluster* vp_cluster = nullptr;
+  if ( true ) { // VP cluster 'hack'
+    auto vp_clusters = getIfExists<LHCb::VPClusters>( LHCb::VPClusterLocation::Default );
+    if ( !vp_clusters || vp_clusters->empty() ) {
+      if ( msgLevel(MSG::DEBUG) ) {
+        debug() << "No VP clusters, skipping..." << endmsg;
+      }
+
+      return StatusCode::SUCCESS;
+    }
+
+    // Find the VP cluster nearest the origin
+    std::vector<std::pair<float,const LHCb::VPCluster*> > scores;
+    for ( const auto& cluster : *vp_clusters ) {
+      scores.emplace_back(
+          std::sqrt( std::pow( cluster->x(), 2) + std::pow( cluster->y(), 2) + std::pow( cluster->z(), 2) ),
+          cluster );
+    }
+
+    auto min_iter = std::min_element( scores.begin(), scores.end() );
+    vp_cluster = min_iter->second;
+  }*/
+
+  const auto& muonTracks = m_trackTool->tracks();
+  auto tracks = std::make_unique<LHCb::Tracks>();
+
+  // -- Loop over all Muon Tracks
+  for ( const auto& muonTrack : muonTracks ) {
+    // -- Get the momentum of the muon track, make a new track
+    auto track = m_momentumTool->recMomentum( const_cast<MuonTrack&>( muonTrack ) ); //@FIXME
+    if ( !track ) {
+      continue;
+    }
+
+    if ( m_fitTracks ) {
+      // Try and improve the covariance information
+      auto sc = m_trackFitter->operator()( *track );
+      if ( sc.isFailure() ) {
+        if ( msgLevel( MSG::WARNING ) ) {
+            warning() << "Track fit failed" << endmsg;
+        }
+
+        continue;
+      }
+    }
+
+    // -- Change q/p until it points to the origin (adapted from Wouter)
+    LHCb::State veloState, muonState;
+    auto sc = iterateToPV(track.get(), muonState, veloState, PVPos, muonTrack.qOverP() );
+    if ( sc.isFailure() ) {
+      Warning( "==> Could not iterate to PV!", StatusCode::SUCCESS, 0 ).ignore();
+      continue;
+    }
+
+    // -- Set Pattern Reco status and track type, finally fit the track
+    track->setPatRecStatus( LHCb::Track::PatRecStatus::PatRecIDs );
+    track->setType( LHCb::Track::Types::Muon );
+
+    if ( !m_fitTracks ) {
+      track->clearStates( ); // remove the state recMomentum created
+    }
+
+    track->addToStates( veloState );
+    track->addToStates( muonState );
+    // ---------------------------------------------------------------
+
+    // -- Finally, insert the track into the container!
+    tracks->insert( track.release() );
+  }
+
+  // ---------------------------------------------------------------
+
+  if(msgLevel(MSG::DEBUG)) {
+    debug() << "Filling " << tracks->size() << " tracks in " << m_outputLoc << endmsg;
+  }
+
+  // Control flow can stop if we didn't make anything
+  setFilterPassed( !tracks->empty() );
+
+  // Put our output on the TES
+  put( tracks.release(), m_outputLoc );
+
+  return StatusCode::SUCCESS;
+}
+
+//=============================================================================
+//  Change the q/p till the track points to the PV (stolen from Wouter)
+//=============================================================================
+StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& muonState, LHCb::State& veloState,
+    const std::vector<double> &PVPos, double qOverP)
+{
+  muonState = track->closestState( 15000 );
+
+  double dXdQoP = 1e7 ;
+  double deltaX = 100;
+  muonState.setQOverP( qOverP ) ;
+
+  // Set the y slope based on the target position at ~the origin
+  muonState.setTy( ( muonState.y() - PVPos[1] ) / ( muonState.z() - PVPos[2] ) );
+
+  // Set the uncertainty on ty to just come from the y uncertainty from the muon stations
+  auto cov = muonState.covariance();
+  cov(3, 3) = std::pow( muonState.ty(), 2 ) * ( cov(1, 1) / ( muonState.y() - PVPos[1] ) );
+  muonState.setCovariance( cov );
+
+  // -- Now call the extrapolator and iterate until we have the desired accuracy.
+  double tolerance = 0.5 ; // [mm]
+  Gaudi::TrackMatrix jacobian ;
+
+  for(int i=0; i<10 && std::abs(deltaX) > tolerance ; ++i) {
+    veloState = muonState;
+    auto sc = m_extrapolator->propagate( veloState, PVPos[2], &jacobian);
+    if ( sc.isFailure() ) {
+      return StatusCode::FAILURE;
+    }
+    dXdQoP = jacobian(0,4) ;
+    deltaX = -(veloState.x() - PVPos[0]) ;
+    double deltaQoP = deltaX / dXdQoP ;
+    muonState.setQOverP( muonState.qOverP() + deltaQoP ) ;
+  }
+
+  veloState.setLocation( LHCb::State::ClosestToBeam );
+  muonState.setLocation( LHCb::State::Muon );
+
+  return m_extrapolator->propagate( veloState, PVPos[2], &jacobian);
+}
diff --git a/Tr/TrackTools/src/MuonSeeding.h b/Tr/TrackTools/src/MuonSeeding.h
new file mode 100644
index 0000000000000000000000000000000000000000..0aafbb802a20d20c8acfb7b2141016a8d0051e4e
--- /dev/null
+++ b/Tr/TrackTools/src/MuonSeeding.h
@@ -0,0 +1,64 @@
+#ifndef MUONSEEDING_H
+#define MUONSEEDING_H 1
+
+// Include files
+// from Gaudi
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "MuonInterfaces/MuonTrack.h"
+#include "MuonInterfaces/MuonHit.h"
+#include "MuonInterfaces/IMuonTrackMomRec.h"
+#include "MuonInterfaces/IMuonTrackRec.h"
+#include "Event/Track.h"
+#include "Event/Particle.h"
+
+#include "Event/State.h"
+#include "Event/RecVertex.h"
+
+#include "TrackInterfaces/ITrackFitter.h"
+#include "TrackInterfaces/ITrackExtrapolator.h"
+
+#include "Kernel/LHCbID.h"
+#include "GaudiKernel/Chrono.h"
+
+
+/** @class MuonSeeding MuonSeeding.h
+ *
+ * \brief  Make a MuonSeeding: Get muon standalone tracks
+ *
+ * Parameters:
+ * - ToolName: Name for the tool that makes muon standalone track.
+ * - Extrapolator: Name for the track extrapolator.
+ * - FillMuonStubInfo: Fill parameters of muon stub in info fields of track;
+ * - Output: The location the tracks should be written to.
+ *
+ *  @author Michel De Cian
+ *  @date   2010-09-20
+ */
+
+class MuonSeeding final : public GaudiAlgorithm {
+public:
+  /// Standard constructor
+  using GaudiAlgorithm::GaudiAlgorithm;
+
+  MuonSeeding( const std::string& name, ISvcLocator* pSvcLocator );
+  StatusCode initialize() override;    ///< Algorithm initialization
+  StatusCode execute() override;    ///< Algorithm execution
+
+private:
+
+  // -- Methods
+  StatusCode fillPVs(std::vector<double>& PVPos);
+  StatusCode iterateToPV(LHCb::Track* track, LHCb::State& muonState, LHCb::State& veloState,
+      const std::vector<double>& PVPos, double qOverP);
+
+  // -- Properties
+  Gaudi::Property<bool>         m_fitTracks         { this, "FitTracks",        true };
+  Gaudi::Property<std::string>  m_outputLoc         { this, "Output",           "Rec/MuonSeeding/Tracks" };
+
+  // -- Tools
+  ToolHandle<IMuonTrackRec>       m_trackTool { "MuonNNetRec/MuonRecTool", this };
+  ToolHandle<ITrackFitter>        m_trackFitter { "TrackMasterFitter/Fitter", this };
+  ToolHandle<IMuonTrackMomRec>    m_momentumTool { "MuonTrackMomRec/MuonMomTool", this };
+  ToolHandle<ITrackExtrapolator>  m_extrapolator { "TrackMasterExtrapolator/Extrapolator", this };
+};
+#endif // GETMUONTRACK_H