diff --git a/Pr/PrAlgorithms/doc/release.notes b/Pr/PrAlgorithms/doc/release.notes
index 8f2dbf7572667853414f8549fe6f65dc6e15d2cd..9b78e544d596bd94c88791f03da089cef6bf53ee 100644
--- a/Pr/PrAlgorithms/doc/release.notes
+++ b/Pr/PrAlgorithms/doc/release.notes
@@ -3,6 +3,10 @@
 ! Responsible : Yasmine Amhis
 ! Purpose     : Pattern reconstruction for the upgrade
 !-----------------------------------------------------------------------------
+
+! 2016-12-14 - Jeroen van Tilburg
+ - Adapt to new FT geometry and numbering scheme
+
 2016-04-18 Renato Quagliani
  - Fix warnings for unused variables
 2016-04-18 Renato Quagliani
diff --git a/Pr/PrAlgorithms/src/PrFTHitManager.cpp b/Pr/PrAlgorithms/src/PrFTHitManager.cpp
index 21d9ab1cabef943f7539cb06db2a67f0decdb9f7..c27bc56980dcf64758c28c38a0a327b9b3f372f7 100644
--- a/Pr/PrAlgorithms/src/PrFTHitManager.cpp
+++ b/Pr/PrAlgorithms/src/PrFTHitManager.cpp
@@ -6,8 +6,6 @@
 // local
 #include "PrFTHitManager.h"
 
-//#include  "Kernel/FTChannelID.h"
-
 //-----------------------------------------------------------------------------
 // Implementation file for class : PrFTHitManager
 //
@@ -17,19 +15,15 @@
 // Declaration of the Tool Factory
 DECLARE_TOOL_FACTORY( PrFTHitManager )
 
-
 //=============================================================================
 // Standard constructor, initializes variables
 //=============================================================================
 PrFTHitManager::PrFTHitManager( const std::string& type,
-                                  const std::string& name,
-                                  const IInterface* parent )
+                                const std::string& name,
+                                const IInterface* parent )
 : PrHitManager ( type, name , parent )
 {
   declareInterface<PrFTHitManager>(this);
-  declareProperty( "DoSmearing" , m_doSmearing = false);
-  declareProperty( "XSmearing",   m_xSmearing  = -1. );
-  declareProperty( "ZSmearing",   m_zSmearing  = -1. );
 }
 
 //=========================================================================
@@ -38,66 +32,43 @@ PrFTHitManager::PrFTHitManager( const std::string& type,
 void PrFTHitManager::buildGeometry ( ) {
   // -- only build geometry once.
   if(m_geometryBuilt) return;
-  if( m_doSmearing){
-    IRndmGenSvc* randSvc = svc<IRndmGenSvc>( "RndmGenSvc", true );
-    StatusCode sc = m_gauss.initialize( randSvc, Rndm::Gauss( 0., 1. ) );
-    if ( sc.isFailure() ){
-      error() << "Could not initialize Rndm::Gauss generator" << endmsg;
-      return;
-    }
-  }
-  m_ftDet = getDet<DeFTDetector>( DeFTDetectorLocation::Default );
 
+  // Retrieve and initialize DeFT (no test: exception in case of failure)
+  m_ftDet = getDet<DeFTDetector>( DeFTDetectorLocation::Default );
   if(UNLIKELY(! m_ftDet ) ){
     error()<<" Cannot load the FTDetector from default location "<<endmsg;
     return;
   }
+  if( m_ftDet->version() < 61 ) {
+    error() << "This version requires FTDet v6.1 or higher" << endmsg;
+    return;
+  }
 
-  //version 20 is the detector with monolayer and fibremat structure, including the dead regions
-  if( msgLevel(MSG::DEBUG) )  debug() << "DETECTOR VERSION: " << m_ftDet->version() << endmsg;
-  //Different Detector Segments has to be used to set the geometry of zone Up and Down to be consistent with new DeFTFibreMat
-  DetectorSegment segUp;
-  DetectorSegment segDown;
-  for ( std::vector<DeFTLayer*>::const_iterator itL = m_ftDet->layers().begin();  //loop over layers
-	m_ftDet->layers().end() != itL; ++itL ) {
-    int id = (*itL)->params()->param<int>("layerID");  //ask layerID
-
-    //FTChannelID (layer,module,mat,SiPMID,netCellID)
-    LHCb::FTChannelID ChUp( id, 2u, 0u, 0u, 0u ); //"0u" means "unsigned 0"
-    const DeFTFibreMat* ftMatUp = m_ftDet->findFibreMat(ChUp);
-    segUp = ftMatUp->createDetSegment( ChUp, 0. );
-
-    LHCb::FTChannelID ChDown( id, 2u, 1u, 0u, 0u ); //"0u" means "unsigned 0"
-    const DeFTFibreMat* ftMatDown = m_ftDet->findFibreMat(ChDown);
-    segDown = ftMatDown->createDetSegment( ChDown, 0. );
-
-    if ( msgLevel(MSG::DEBUG) ){
-      debug() << "STEREO ANGLE Up: " <<  ftMatUp->angle() << endmsg;
-      debug() << "STEREO ANGLE Down: " <<  ftMatDown->angle() << endmsg;
-    }
-
-
-    //The setGeometry defines the z at y=0, the dxDy and the dzDy, as well as the isX properties of the zone.
-    //This is important, since these are used in the following.
-    //They are set once for each zone in this method.
-    MakeZone( 2*id, segUp , -4090., 4090., -50., 3030. );
-    MakeZone( 2*id+1 , segDown, -4090., 4090., -3030., 50. );
-
-    //----> Debug zones
-    if( msgLevel(MSG::DEBUG)){
-      debug() << "Layer " << id << " z " << zone(2*id).z()
-              << " angle " << zone(2*id).dxDy() << endmsg;
+  // Loop over all layers
+  for( auto station : m_ftDet->stations() ) {
+    for( auto layer : station->layers() ) {
+      int id = 4*(station->stationID() - 1) + layer->layerID();
+
+      DetectorSegment seg(0, layer->globalZ(), layer->dxdy(), layer->dzdy(), 0., 0.);
+      float xmax = 0.5*layer->sizeX();
+      float ymax = 0.5*layer->sizeY();
+
+      //The setGeometry defines the z at y=0, the dxDy and the dzDy, as well as the isX properties of the zone.
+      //This is important, since these are used in the following.
+      //They are set once for each zone in this method.
+      MakeZone( 2*id,   seg, -xmax, xmax, -25., ymax ); // Small overlap (25 mm) for stereo layers
+      MakeZone( 2*id+1, seg, -xmax, xmax, -ymax, 25. ); // Small overlap (25 mm) for stereo layers
+
+      //----> Debug zones
+      if( msgLevel(MSG::DEBUG)){
+        debug() << "Layer " << id << " z " << zone(2*id).z()
+                << " angle " << zone(2*id).dxDy() << endmsg;
+      }
+      //----> Debug zones
     }
-    //----> Debug zones
   }
-
-  //---> Debug smearing
-  if( msgLevel(MSG::DEBUG)) debug() << "XSmearing " << m_xSmearing << " ZSmearing " << m_zSmearing << endmsg;
-  //---> Debug smearing
-
+  
   m_geometryBuilt = true;
-
-
 }
 
 
@@ -107,63 +78,41 @@ void PrFTHitManager::buildGeometry ( ) {
 void PrFTHitManager::decodeData ( ) {
 
   if(m_decodedData) return;
-
-  if( msgLevel(MSG::DEBUG) ) debug() << "I AM IN DECODEDATA " << endmsg;
-
+  
   typedef FastClusterContainer<LHCb::FTLiteCluster,int> FTLiteClusters;
-  FTLiteClusters* clus = get<FTLiteClusters>( LHCb::FTLiteClusterLocation::Default );
-  if( msgLevel(MSG::DEBUG) ) debug() << "Retrieved " << clus->size() << " clusters" << endmsg;
-  const DeFTFibreMat* ftMat = nullptr;
-  const DeFTFibreMat* anaFtMat = nullptr;
-  unsigned int oldFibreMatID = 99999999;
-
-  for ( FTLiteClusters::iterator itC = clus->begin(); clus->end() != itC; ++itC ) {
-    /// find fibremat to which the cluster belongs
-    anaFtMat = m_ftDet->findFibreMat( (*itC).channelID() );
-    if(anaFtMat->FibreMatID() != oldFibreMatID)  {
-      oldFibreMatID =  anaFtMat->FibreMatID();
-      ftMat =  anaFtMat;
-      if ( nullptr == ftMat ) {
-        info() << "FiberMat not found for FT channelID " << (*itC).channelID() << endmsg;
-        oldFibreMatID = 99999999;
-      }
+  FTLiteClusters* clusters = get<FTLiteClusters>( LHCb::FTLiteClusterLocation::Default );
+  if ( msgLevel( MSG::DEBUG) ) debug() << "Retrieved " << clusters->size() << " clusters" << endmsg;
+
+  const DeFTMat* mat = nullptr;
+  unsigned int prevMatID = 99999999;
+  for ( auto clus : *clusters ) {
+    if( clus.channelID().uniqueMat() != prevMatID ) {
+      mat = m_ftDet->findMat( clus.channelID() );
+      prevMatID = mat->elementID().uniqueMat();
     }
 
-    double fraction = (*itC).fraction() + 0.125;   // Truncated to 4 bits = 0.25. Add half of it
-    LHCb::FTChannelID id = (*itC).channelID();
-    int lay  = (*itC).channelID().layer();
-    int zone = (int)((*itC).channelID().mat() > 0);  //if yes, it is bottom (now bottom is 1)
+    double fraction = clus.fraction();
+    LHCb::FTChannelID id = clus.channelID();
+
+    int lay  = 4*(id.station()-1) + id.layer();
+    int zone = int(mat->isBottom());
     int code = 2*lay + zone;
 
-    //variables
-    const double fibreMatHalfSizeY = ftMat->getFibreMatHalfSizeY();
-    const double cellSizeX = ftMat->getCellSizeX();
-    const float dzDy = ftMat->getDzDy();
-    const float dxDy = -ftMat->getTanAngle();// = -m_tanAngle
-    const double hitLocalX = ftMat->cellLocalX(id) + fraction*cellSizeX;
-    double xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz;
-    ftMat->geometry()->toGlobalMatrix().GetComponents( xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz );
-
-    //const float yBorder = yx* hitLocalX + yy * (2*zone-1)*fibreMatHalfSizeY + dy;
-    const float yMin = yx* hitLocalX - yy *fibreMatHalfSizeY + dy;
-    const float yMax = yx* hitLocalX + yy *fibreMatHalfSizeY + dy;
-
-    //calculate hit x,z at y=0 in global frame (as local y=0 NOT global y=0)
-    const double hitGlobala_X = hitLocalX * xx + dx;
-    const double hitGlobalb_X = hitGlobala_X + fibreMatHalfSizeY*xy;
-    const double hitGlobala_Y = hitLocalX * yx+ dy;
-    const double hitGlobalb_Y = hitGlobala_Y + fibreMatHalfSizeY*yy;
-    const double hitGlobala_Z = hitLocalX * zx + dz;
-    const double hitGlobalb_Z = hitGlobala_Z + fibreMatHalfSizeY*zy;
-    const double axy = (hitGlobala_X-hitGlobalb_X)/(hitGlobala_Y-hitGlobalb_Y);
-    const double azy = (hitGlobala_Z-hitGlobalb_Z)/(hitGlobala_Y-hitGlobalb_Y);
-    const float x0 = hitGlobalb_X-axy*hitGlobalb_Y;
-    const float z0 = hitGlobalb_Z-azy*hitGlobalb_Y;
-    //Error X position
-    float errX = 0.05 + .03 * (*itC).size();
-    addHitInZone( code, LHCb::LHCbID( (*itC).channelID() ), x0, z0, dxDy, dzDy, yMin, yMax, errX , zone, lay );
-    if ( msgLevel(MSG::DEBUG) ) info() << " .. hit " << (*itC).channelID() << " zone " << zone << " x " << x0 << endmsg;
+    auto endPoints = mat->endPoints(id, fraction);
+    double invdy = 1./(endPoints.first.y()-endPoints.second.y());
+    double dxdy  = (endPoints.first.x()-endPoints.second.x())*invdy;
+    double dzdy  = (endPoints.first.z()-endPoints.second.z())*invdy;
+    float x0     = endPoints.first.x()-dxdy*endPoints.first.y();
+    float z0     = endPoints.first.z()-dzdy*endPoints.first.y();
+    float yMin   = std::min(endPoints.first.y(), endPoints.second.y());
+    float yMax   = std::max(endPoints.first.y(), endPoints.second.y());
+
+    float errX = 0.170; // TODO: this should be ~80 micron; get this from a tool
+    addHitInZone( code, LHCb::LHCbID( id ), x0, z0, dxdy, dzdy, yMin, yMax, errX , zone, lay );
+    if ( msgLevel(MSG::VERBOSE) )
+      verbose() << " .. hit " << id << " code=" << code << " x=" << x0 << " z=" << z0 << endmsg;
   }
+
   for ( unsigned int lay = 0; nbZones() > lay ; ++lay ) {
     std::sort( getIterator_Begin(lay), getIterator_End(lay), PrHit::LowerByX0() );
   }
diff --git a/Pr/PrAlgorithms/src/PrFTHitManager.h b/Pr/PrAlgorithms/src/PrFTHitManager.h
index 8e9ad8b1021342d8779842333bd54874d0d1c62a..7a84e409b0bede5ad13e3ee98c5e7cf12dc8d038 100644
--- a/Pr/PrAlgorithms/src/PrFTHitManager.h
+++ b/Pr/PrAlgorithms/src/PrFTHitManager.h
@@ -13,10 +13,6 @@ static const InterfaceID IID_PrFTHitManager ( "PrFTHitManager", 1, 0 );
  *  Tool that transforms clusters into 'hits' (spatial positions) which are then used by the pattern
  *  recognition algorithms involving the FT.
  *
- *  Parameters:
- *  - XSmearing: Amount of gaussian smearing in x
- *  - ZSmearing: Switch on displacement in z (not used at the moment)
- *
  *  @author Olivier Callot
  *  @date   2012-03-13
  */
@@ -38,14 +34,9 @@ public:
    */
   void decodeData() override;
 
-
 protected:
 
  private:
-  DeFTDetector* m_ftDet;
-  bool m_doSmearing = false;
-  Rndm::Numbers m_gauss; ///< Random number generator for gaussian smearing
-  float m_xSmearing =-1.;     ///< Amount of smearing in x
-  float m_zSmearing =-1.; 
+  DeFTDetector* m_ftDet = nullptr;
 };
 #endif // PRFTHITMANAGER_H
diff --git a/Pr/PrMCTools/CMakeLists.txt b/Pr/PrMCTools/CMakeLists.txt
index 82046778da0c73b4ae67e26e141ca32a76e53161..a682e63eeef664f02075552780a0e383abee2d1f 100644
--- a/Pr/PrMCTools/CMakeLists.txt
+++ b/Pr/PrMCTools/CMakeLists.txt
@@ -3,7 +3,8 @@
 ################################################################################
 gaudi_subdir(PrMCTools v2r12)
 
-gaudi_depends_on_subdirs(Det/VPDet
+gaudi_depends_on_subdirs(Det/FTDet       
+                         Det/VPDet
                          Det/STDet
                          Event/FTEvent
                          Event/LinkerEvent
@@ -19,5 +20,5 @@ include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS})
 gaudi_add_module(PrMCTools
                  src/*.cpp
                  INCLUDE_DIRS Event/FTEvent Pr/PrPixel
-                 LINK_LIBRARIES LinkerEvent MCEvent GaudiAlgLib PrKernel VPDetLib STDetLib LoKiMCLib)
+                 LINK_LIBRARIES LinkerEvent MCEvent GaudiAlgLib PrKernel VPDetLib STDetLib FTDetLib LoKiMCLib)
 
diff --git a/Pr/PrMCTools/doc/release.notes b/Pr/PrMCTools/doc/release.notes
index 93dd0ad667a418f32e0318397a01519152855a37..429d61f8032e048c22592422af169d59e162264e 100644
--- a/Pr/PrMCTools/doc/release.notes
+++ b/Pr/PrMCTools/doc/release.notes
@@ -4,6 +4,10 @@
 ! Purpose     : MC tools for pattern recognition for the Upgrade
 !-----------------------------------------------------------------------------
 
+! 2016-12-14 - Jeroen van Tilburg
+ - Adapt to new FT geometry and numbering scheme
+ - Changed location of Linker tables. 
+
 !========================= PrMCTools v2r12 2016-01-28 =========================
 ! 2016-01-21 Marco Cattaneo
  - Use std::abs instead of fabs for int argument, to silence clang warning
diff --git a/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp b/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp
index 061c6a2bd81460305487cf02ec5da0bc4f1624d1..7d1d9575b79d905ad0695e19df4041575020c326 100644
--- a/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp
+++ b/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp
@@ -4,7 +4,7 @@
 #include "Event/Track.h"
 #include "Event/StateParameters.h"
 #include "PrCheatedSciFiTracking.h"
-#include "Event/FTCluster.h"
+#include "Event/FTLiteCluster.h"
 #include "Linker/LinkedTo.h"
 #include "Event/MCTrackInfo.h"
 
@@ -93,7 +93,7 @@ StatusCode PrCheatedSciFiTracking::finalize() {
 //=========================================================================
 void PrCheatedSciFiTracking::makeLHCbTracks ( LHCb::Tracks* result ) {
 
-  LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
+  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
   LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default );
 
   MCTrackInfo trackInfo( evtSvc(), msgSvc() );
diff --git a/Pr/PrMCTools/src/PrClustersResidual.cpp b/Pr/PrMCTools/src/PrClustersResidual.cpp
index 56c3a23fcd15a828fa3b5d527ed5b1de0668102a..a1b5fc54bc7f401d47161f2c95e9fbbf909d9e94 100644
--- a/Pr/PrMCTools/src/PrClustersResidual.cpp
+++ b/Pr/PrMCTools/src/PrClustersResidual.cpp
@@ -38,16 +38,6 @@ PrClustersResidual::PrClustersResidual( const std::string& name,
 : GaudiTupleAlg ( name , pSvcLocator ),
   m_ftHitManager(nullptr),
   m_zone(24){
-  declareProperty("MCHitsLocation",m_mcHitLocation = "/Event/MC/FT/Hits");
-  declareProperty("HitManagerName",m_hitManagerName = "PrFTHitManager");
-  declareProperty("DebugTracking",m_debugTracking = false);
-  declareProperty("DoClusterResidual",m_doClusterResidual = false); //Cluster residual tuple
-  declareProperty("DoTrackStudy",m_doTrackStudy = false); // Produce tuple for track studies
-  declareProperty("DecodeData",m_decodeData = false); //ask to decode data ( is False if other algorithms are runned
-  declareProperty("OnlyHasT",m_onlyHasT=false); //(keep only MCParticles with hasT = true
-  declareProperty("RemoveClones",m_removeClones = true); //Observed clones MCHits when processed ( not understood )
-  declareProperty("DumpAllHits",m_dumpAllHits = false);
-  declareProperty("Occupancy",m_Occupancy = true);
 }
 //=============================================================================
 // Destructor
@@ -83,8 +73,17 @@ StatusCode PrClustersResidual::initialize() {
 //=============================================================================
 // Main execution
 //=============================================================================
-StatusCode PrClustersResidual::execute() {
-
+StatusCode PrClustersResidual::execute() {  
+  m_ftDet = getDet<DeFTDetector>(DeFTDetectorLocation::Default);
+  if( UNLIKELY( !m_ftDet) ){
+    error() << " Cannot load the FTDetector from default location" << endmsg;
+    return StatusCode::FAILURE;
+  }else if( m_ftDet->version() < 61){
+    error()<< "This code requires FTDet v6.1 or higher "<<endmsg;
+    return StatusCode::FAILURE;
+  }
+  
+  
   if ( msgLevel(MSG::DEBUG) ) debug() << "==> Execute" << endmsg;
   m_nEvents++;
   if(m_doClusterResidual){
@@ -124,8 +123,8 @@ void PrClustersResidual::Occupancy(){
   //"/Event/Prev/MC/FT/Hits"
 
   LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" );
-  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(),   LHCb::FTClusterLocation::Default );
-  LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits");
+  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), m_ftClusterToParticleLinkerLocation );
+  LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), m_ftClusterToHitLinkerLocation);
   char layerName[100];
   char Title[100];
   std::vector<int> nHits(12,0);
@@ -149,7 +148,7 @@ void PrClustersResidual::Occupancy(){
   int recoblewanted = 0;
   LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default );
   LHCb::MCParticle* mcPart = nullptr;
-  for( LHCb::MCParticles::const_iterator iPart = mcParts->begin(); iPart != mcParts->end(); ++iPart){
+  for( auto iPart = mcParts->begin(); iPart != mcParts->end(); ++iPart){
     mcPart = (*iPart);
     bool isRecoble = trackInfo.hasT(mcPart);
     bool recobleWanted = isRecoble && (mcPart->originVertex()->isPrimary() || mcPart->originVertex()->isDecay());
@@ -174,19 +173,30 @@ void PrClustersResidual::Occupancy(){
     for(auto itH = m_ftHitManager->getIterator_Begin(zoneI); m_ftHitManager->getIterator_End(zoneI) != itH;++itH){
       hit = &(*itH);
       nHits[layer]++;
-      LHCb::FTLiteCluster liteCluster = getLiteCluster(hit->id());
+      auto liteCluster = getLiteCluster(hit->id());
+      
+      //---- new numbering scheme (v61 sci-fi) plots for occupancy! pseudochannel conversion needed!
+      /*
+        It requires to have the m_ftDet as private member of the algorithm and also modify the CMakeLists to compile
+        PrMCTools linking the Det/FTDet package.
+      */
+      int quarter = liteCluster.channelID().quarter();
+      const auto chanID = liteCluster.channelID();
+      const DeFTModule* module = m_ftDet->findModule( chanID );
+      if( module == nullptr){
+        error()<<"Cannot find Module associated to FTChannelID : "<<chanID<<endmsg;
+        error()<<"Abort Occupancy plot checking"<<endmsg;
+        return;
+        }
+      int pseudochannel = module->pseudoChannel( chanID );
+      //--- runs from [0 to 12,287] (96 SiPm *128 channels/SiPm)
+      //--- 4 quartes [0,1,2,3] , 128 channels squeezed to 1 value. (Occupancy in terms of Hits/SiPM array [128 channels]
+      //--- ToDo : make it depenent on xml geometry information through m_ftDet 
+      int sipm = std::floor( pseudochannel/128);
+      int index = quarter * 128  +  sipm ;
+      // old way [ < v61 ] : int index = liteCluster.channelID().uniqueSiPM() & 511;
+      
       
-      int Module = liteCluster.channelID().module();
-      int Quarter = liteCluster.channelID().quarter();
-      int SipmID = liteCluster.channelID().sipmId();
-      int layer = liteCluster.channelID().layer();
-      int module2 = Module;
-      if(Module == 11 && Module > 4 ) module2 = 1;
-      if(Module != 11 && Module > 4 ) module2 = Module-3;
-      if(Module <  5  ) module2 = 6-Module;
-      if(Module == 10 ) module2 = 1;
-      //---- VALID UNTIL v5x9 SciFi geometry. From v6x1, change in numbering scheme leads to different indexing.
-      int index = Quarter*128 + 16*module2+ std::floor(SipmID);
       sprintf(layerName,"Layer%i/AllContributions",layer);
       sprintf(Title,"Sipm Index Layer %i;SiPm Id;Clusters/Sipm/Event",layer);
       plot1D(index,layerName,Title,-0.5,511.5,512);
@@ -308,8 +318,8 @@ void PrClustersResidual::TrackStudy(){
   LHCb::ODIN* odin=getIfExists<LHCb::ODIN>(LHCb::ODINLocation::Default);
   //Linker for the MCHit
   LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" );
-  LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
-  LinkedTo<LHCb::MCHit, LHCb::FTCluster> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits");
+  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
+  LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits");
   //Get the linker for the generation of the track to study (the one after the linker)
   //Get the McParticles
   LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default );
@@ -430,6 +440,8 @@ void PrClustersResidual::TrackStudy(){
       debug()<<"For the processed MCParticle i found  "<<track.size()<<"PrHits   associated to "<<assoc_MCHit.size()<<"   MCHits"<<endmsg;
       tuple->farray("FiredLayers_Counter",firedLayers,"FiredLayers",24);
       tuple->column("CheatedSeeding_NHits",totHits);
+      std::vector<double> MCHit_EntryX,MCHit_EntryY,MCHit_EntryZ;      
+      std::vector<double> MCHit_ExitX,MCHit_ExitY,MCHit_ExitZ;
       std::vector<double> MCHit_X,MCHit_Y,MCHit_Z;
       std::vector<double> MCHit_tx,MCHit_ty,MCHit_p;
       std::vector<double> MCHit_pathlength;
@@ -443,12 +455,12 @@ void PrClustersResidual::TrackStudy(){
 
       //ChannelID of track hits
       std::vector<double> ChID_Fraction;
-      std::vector<double> ChID_Charge;
       std::vector<double> ChID_SipmCell;
+      std::vector<double> ChID_Mat;
       std::vector<double> ChID_Module;
       std::vector<double> ChID_Layer;
       std::vector<double> ChID_Quarter;
-      std::vector<double> ChID_Mat;
+      std::vector<double> ChID_Station;
       std::vector<double> ChID_SipmID;
       std::vector<double> ChID_ID;
 
@@ -477,6 +489,13 @@ void PrClustersResidual::TrackStudy(){
           MCHit_ty.push_back(mcHit->dydz());
           MCHit_p.push_back(mcHit->p());
           MCHit_pathlength.push_back(mcHit->pathLength());
+          
+          MCHit_EntryX.push_back(mcHit->entry().X());
+          MCHit_EntryY.push_back(mcHit->entry().Y());
+          MCHit_EntryZ.push_back(mcHit->entry().Z());
+          MCHit_ExitX.push_back(mcHit->exit().X());
+          MCHit_ExitY.push_back(mcHit->exit().Y());
+          MCHit_ExitZ.push_back(mcHit->exit().Z());
           MCHit_X.push_back(mcHit->midPoint().X());
           MCHit_Y.push_back(mcHit->midPoint().Y());
           MCHit_Z.push_back(mcHit->midPoint().Z());
@@ -501,14 +520,14 @@ void PrClustersResidual::TrackStudy(){
             }
             LHCb::FTLiteCluster liteCluster = getLiteCluster(hit->id());
             ChID_Fraction.push_back(liteCluster.fraction());
-            ChID_SipmCell.push_back(liteCluster.channelID().sipmCell());
+            ChID_SipmCell.push_back(liteCluster.channelID().channel());
             ChID_Module.push_back(liteCluster.channelID().module());
+            ChID_Mat.push_back(liteCluster.channelID().mat());
             ChID_Layer.push_back(liteCluster.channelID().layer());
             ChID_Quarter.push_back(liteCluster.channelID().quarter());
-            ChID_Mat.push_back(liteCluster.channelID().mat());
-            ChID_SipmID.push_back(liteCluster.channelID().sipmId());
+            ChID_Station.push_back(liteCluster.channelID().station());
+            ChID_SipmID.push_back(liteCluster.channelID().sipm());
             ChID_ID.push_back(liteCluster.channelID());
-            ChID_Charge.push_back(liteCluster.charge());
             PrHit_LHCBID.push_back(hit->id().ftID());
             PrHit_Xat0.push_back(hit->x());
             PrHit_Zat0.push_back(hit->z());
@@ -530,6 +549,12 @@ void PrClustersResidual::TrackStudy(){
         MCHit_time.push_back(-999999.);
         MCHit_Particle_key.push_back(-999999.);
         MCHit_pathlength.push_back(-999999.);
+        MCHit_EntryX.push_back(-999999.);
+        MCHit_EntryY.push_back(-999999.);
+        MCHit_EntryZ.push_back(-999999.);
+        MCHit_ExitX.push_back(-999999.);
+        MCHit_ExitY.push_back(-999999.);
+        MCHit_ExitZ.push_back(-999999.);
         MCHit_X.push_back(-999999.);
         MCHit_Y.push_back(-999999.);
         MCHit_Z.push_back(-999999.);
@@ -552,6 +577,13 @@ void PrClustersResidual::TrackStudy(){
       tuple->farray("MCHit_tx",MCHit_tx,"MC_ass",100);
       tuple->farray("MCHit_p",MCHit_p,"MC_ass",100);
       tuple->farray("MCHit_pathlength",MCHit_pathlength,"MC_ass",100);
+      tuple->farray("MCHit_Assoc_EntryX",MCHit_EntryX ,"MC_ass",100);
+      tuple->farray("MCHit_Assoc_EntryY",MCHit_EntryY ,"MC_ass",100);
+      tuple->farray("MCHit_Assoc_EntryZ",MCHit_EntryZ ,"MC_ass",100);
+      tuple->farray("MCHit_Assoc_ExitX",MCHit_ExitX ,"MC_ass",100);
+      tuple->farray("MCHit_Assoc_ExitY",MCHit_ExitY ,"MC_ass",100);
+      tuple->farray("MCHit_Assoc_ExitZ",MCHit_ExitZ ,"MC_ass",100);
+
       tuple->farray("MCHit_Assoc_X",MCHit_X ,"MC_ass",100);
       tuple->farray("MCHit_Assoc_Y",MCHit_Y ,"MC_ass",100);
       tuple->farray("MCHit_Assoc_Z",MCHit_Z ,"MC_ass",100);
@@ -571,13 +603,13 @@ void PrClustersResidual::TrackStudy(){
       tuple->farray("PrHit_dzDy",PrHit_dzDy,"PrHit",100);
 
       tuple->farray("ChID_Fraction",ChID_Fraction,"ChID",100);
-      tuple->farray("ChID_SipmCell",ChID_SipmCell,"ChID",100);
-      tuple->farray("ChID_Module",ChID_Module,"ChID",100);
-      tuple->farray("ChID_Layer",ChID_Layer,"ChID",100);
-      tuple->farray("ChID_Quarter",ChID_Quarter,"ChID",100);
-      tuple->farray("ChID_Charge",ChID_Charge,"ChID",100);
-      tuple->farray("ChID_Mat",ChID_Mat,"ChID",100);
+      tuple->farray("ChID_SipmCell",ChID_SipmCell,"ChID",100); 
       tuple->farray("ChID_SipmID",ChID_SipmID,"ChID",100);
+      tuple->farray("ChID_Mat",ChID_Mat,"ChID",100);
+      tuple->farray("ChID_Module",ChID_Module,"ChID",100); 
+      tuple->farray("ChID_Quarter",ChID_Quarter,"ChID",100);
+      tuple->farray("ChID_Layer",ChID_Layer,"ChID",100); 
+      tuple->farray("ChID_Station",ChID_Station,"ChID",100);
       tuple->farray("ChID_ID",ChID_ID,"ChID",100);
 
       tuple->column("N_PrHit_Assoc",track.size());
@@ -958,10 +990,10 @@ void  PrClustersResidual::ClusterResidual(){
   //LinkedTo<LHCb::MCParticle, LHCb::Track> myForwardLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Forward);
   //LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
   if(msgLevel(MSG::DEBUG)) debug()<<" Get ClusterLinker "<<endmsg;
-  LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
+  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
   if(msgLevel(MSG::DEBUG)) debug()<<" Get Cluster to MCHit Linker "<<endmsg;
-
-  LinkedTo<LHCb::MCHit, LHCb::FTCluster> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits");
+                  
+  LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits");
   MCTrackInfo trackInfo( evtSvc(), msgSvc() );
   if(msgLevel(MSG::DEBUG)) debug()<<" Add MCHitAndTrackStudy "<<endmsg;
   Tuples::Tuple tuple = GaudiTupleAlg::nTuple("ClusterMCHitAndTrackStudy","Events");
@@ -1032,6 +1064,8 @@ void  PrClustersResidual::ClusterResidual(){
       if(msgLevel(MSG::DEBUG)) debug()<<"Loaded"<<endmsg;
       //took first MCHit associated to Hit channelID
       std::vector<double>  MCHit_X,MCHit_Y,MCHit_Z;
+      std::vector<double>  MCHit_EntryX,  MCHit_EntryY, MCHit_EntryZ;
+      std::vector<double>  MCHit_ExitX,   MCHit_ExitY,  MCHit_ExitZ;
       std::vector<double>  Residual_X,Residual_Z;
       std::vector<double>  SlopeX;
       std::vector<double>  SlopeY;
@@ -1047,15 +1081,25 @@ void  PrClustersResidual::ClusterResidual(){
         MCHit_X.push_back(-9999.);
         MCHit_Y.push_back(-9999.);
         MCHit_Z.push_back(-9999.);
-        SlopeX.push_back(-999.);
-        SlopeY.push_back(-999.);
-        pathlength.push_back(-999.);
-        energy.push_back( -999.);
+        MCHit_EntryX.push_back(-9999.);
+        MCHit_EntryY.push_back(-9999.);
+        MCHit_EntryZ.push_back(-9999.);
+        
+        MCHit_ExitX.push_back(-9999.);
+        MCHit_ExitY.push_back(-9999.);
+        MCHit_ExitZ.push_back(-9999.);
+
+        SlopeX.push_back(-9999.);
+        SlopeY.push_back(-9999.);
+        pathlength.push_back(-9999.);
+        energy.push_back( -9999.);
         hit_Physics.push_back(false);
       }
       while(mcHit != nullptr){
         nMCHitToCluster++;
-        Gaudi::XYZPoint pMid = mcHit->midPoint();
+        auto pMid = mcHit->midPoint();
+        auto pEntry = mcHit->entry();
+        auto pExit  = mcHit->exit();
         Residual_X.push_back( (*iHit).x(pMid.y()) - pMid.x()); //x at y - y at z
         Residual_Z.push_back( (*iHit).z(pMid.y()) - pMid.z());
         SlopeX.push_back( mcHit->dxdz());
@@ -1065,6 +1109,13 @@ void  PrClustersResidual::ClusterResidual(){
         MCHit_X.push_back( pMid.x());
         MCHit_Y.push_back( pMid.y());
         MCHit_Z.push_back( pMid.z());
+        
+        MCHit_ExitX.push_back( pExit.x());
+        MCHit_ExitY.push_back( pExit.y());
+        MCHit_ExitZ.push_back( pExit.z());
+        MCHit_EntryX.push_back( pEntry.x());
+        MCHit_EntryY.push_back( pEntry.y());
+        MCHit_EntryZ.push_back( pEntry.z());
         //mcHit = myFTCluster2MCHitLink.next();
         // const LHCb::MCParticle *mcPart = mcHit->mcParticle();
         const LHCb::MCParticle *mcPart = mcHit->mcParticle();
@@ -1074,6 +1125,7 @@ void  PrClustersResidual::ClusterResidual(){
         if(mcPart!= nullptr){
           if(msgLevel(MSG::DEBUG)) debug() <<"Getting Vertex"<<endmsg;
           const LHCb::MCVertex* vertex = mcPart->originVertex();
+          //track of physical interest if the associated vertex is a Primary one or from Decay ( i.e. veto the interaction one )
           if( vertex!=nullptr && (vertex->isPrimary() || vertex->isDecay())){
             hit_Physics.push_back( true );
           }else{
@@ -1084,7 +1136,6 @@ void  PrClustersResidual::ClusterResidual(){
       }
       if(msgLevel(MSG::DEBUG)) debug()<<"Storing vectors"<<endmsg;
       tuple->column("isNoiseCluster",isNoise);
-
       tuple->farray("MCHit_Physics", hit_Physics.begin(),hit_Physics.end(),"N_MC",100);
       tuple->farray("CosSlopeX",SlopeX.begin(),SlopeX.end(),"N_MC",100);
       tuple->farray("CosSlopeY",SlopeY.begin(),SlopeY.end(),"N_MC",100);
@@ -1095,17 +1146,26 @@ void  PrClustersResidual::ClusterResidual(){
       tuple->farray("MCHit_X",MCHit_X.begin(),MCHit_X.end(),"N_MC",100);
       tuple->farray("MCHit_Z",MCHit_Z.begin(),MCHit_Z.end(),"N_MC",100);
       tuple->farray("MCHit_Y",MCHit_Y.begin(),MCHit_Y.end(),"N_MC",100);
+      
+      tuple->farray("MCHit_EntryX",MCHit_EntryX.begin(),MCHit_EntryX.end(),"N_MC",100);
+      tuple->farray("MCHit_EntryZ",MCHit_EntryZ.begin(),MCHit_EntryZ.end(),"N_MC",100);
+      tuple->farray("MCHit_EntryY",MCHit_EntryY.begin(),MCHit_EntryY.end(),"N_MC",100);
+      tuple->farray("MCHit_ExitX",MCHit_ExitX.begin(),MCHit_ExitX.end(),"N_MC",100);
+      tuple->farray("MCHit_ExitZ",MCHit_ExitZ.begin(),MCHit_ExitZ.end(),"N_MC",100);
+      tuple->farray("MCHit_ExitY",MCHit_ExitY.begin(),MCHit_ExitY.end(),"N_MC",100);
+
+
       //here fill the tuple with watever you want
       //get Cluster associated to the Hit
       if(msgLevel(MSG::DEBUG)) debug()<<"Getting LiteCluster from PrHit"<<endmsg;
       LHCb::FTLiteCluster litecluster = getLiteCluster( (*iHit).id());
-      tuple->column("ClusterCharge",litecluster.charge());
-      tuple->column("ClusterSize",litecluster.size());
+      tuple->column("ClusterSize",litecluster.isLarge());
       tuple->column("ClusterFraction",litecluster.fraction());
       tuple->column("ClusterChannelID",litecluster.channelID());
-      tuple->column("ClusterChannelIDSipmCell",litecluster.channelID().sipmCell());
-      tuple->column("ClusterChannelIDSipmID",litecluster.channelID().sipmId());
+      tuple->column("ClusterChannelIDSipmCell",litecluster.channelID().channel());
+      tuple->column("ClusterChannelIDSipmID",litecluster.channelID().sipm());
       tuple->column("ClusterChannelIDMat",litecluster.channelID().mat());
+      tuple->column("ClusterChannelIDStation",litecluster.channelID().station());
       tuple->column("ClusterChannelIDModule",litecluster.channelID().module());
       tuple->column("ClusterChannelLayer",litecluster.channelID().layer());
       tuple->column("ClusterChannelQuarter",litecluster.channelID().quarter());
@@ -1152,8 +1212,8 @@ void  PrClustersResidual::ClusterResidual(){
       tuple->column("Hit_Zone",(*iHit).zone());
       tuple->column("Hit_dzDy_manually",((*iHit).z(1000.)-(*iHit).z(0.))/1000);
       //Get the list of tracks generated by the Hit i am lookign to
-      if(m_debugTracking)
-      {
+      if(m_debugTracking){
+        //get linker Forward tracks and Seeding ones to see hit content!
         LinkedTo<LHCb::MCParticle, LHCb::Track> myForwardLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Forward);
         LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
 
@@ -1311,3 +1371,4 @@ LHCb::FTLiteCluster PrClustersResidual::getLiteCluster(const LHCb::LHCbID id)
 
   return idTracks;
 }
+
diff --git a/Pr/PrMCTools/src/PrClustersResidual.h b/Pr/PrMCTools/src/PrClustersResidual.h
index 0c6c12778d4656d1cc176c85f6927e87be051a29..6c48866a672273afd88ab593cacfb9d8dda1680a 100644
--- a/Pr/PrMCTools/src/PrClustersResidual.h
+++ b/Pr/PrMCTools/src/PrClustersResidual.h
@@ -23,7 +23,7 @@
 #include "Event/Track.h"
 #include "Linker/LinkerTool.h"
 #include "PrKernel/IPrCounter.h"
-
+#include "FTDet/DeFTDetector.h"
 /** @class PrClustersResidual PrClustersResidual.h
  *  Make nTuples of Cluster vs MCHits and Hit content on track + Occupancy study possible ( also in PrPlotFTHits )
  *  -MCHitsLocation : Where to find the MCHits ( /Event/MC/FT/FTHits by default )
@@ -103,37 +103,28 @@ private:
   std::vector<const LHCb::Track*> getTrack(const LHCb::LHCbID id, const std::string location);
   LHCb::FTLiteCluster getLiteCluster(const LHCb::LHCbID id);
   
-  std::string m_mcHitLocation;
   
-  //  unsigned int m_zone;
-  bool m_debugTracking;
+  Gaudi::Property<std::string> m_mcHitLocation                     {this, "MCHitsLocation", "/Event/MC/FT/Hits"};
+  Gaudi::Property<std::string> m_ftClusterToParticleLinkerLocation { this, "FTClusterToParticleLinkerLocation", LHCb::FTLiteClusterLocation::Default + "WithSpillover"};
+  Gaudi::Property<std::string> m_ftClusterToHitLinkerLocation      { this, "FTClusterToHitLinkerLocation", LHCb::FTLiteClusterLocation::Default + "2MCHitsWithSpillover"};
+  Gaudi::Property<std::string> m_hitManagerName                     { this, "HitManagerName", "PrFTHitManager"};
   
+  Gaudi::Property<bool> m_debugTracking                             { this, "DebugTracking", false};
+  Gaudi::Property<bool> m_doClusterResidual                         { this, "DoClusterResidual", false}; //Cluster residual tuple
+  Gaudi::Property<bool> m_doTrackStudy                              { this, "DoTrackStudy", false}; // Produce tuple for track studies
+  Gaudi::Property<bool> m_decodeData                                { this, "DecodeData", false}; //ask to decode data ( is False if other algorithms are runned
+  Gaudi::Property<bool> m_onlyHasT                                  { this, "OnlyHasT",false}; //(keep only MCParticles with hasT = true
+  Gaudi::Property<bool> m_removeClones                              { this, "RemoveClones", true}; //Observed clones MCHits when processed ( not understood )
+  Gaudi::Property<bool> m_dumpAllHits                               { this, "DumpAllHits", false }; //Dump in tuple for a given event all the hits
+  Gaudi::Property<bool> m_Occupancy                                 { this, "Occupancy",   true };// Produce Occupancy plot
   int m_nEvents;
-  
   double m_NClusters ;
   double m_NMCHit ;
   double m_NMCHit_inClusters;
-  bool m_onlyHasT;
   PrHitManager* m_ftHitManager;
   unsigned int m_zone;
   
-  bool m_decodeData;
-  bool m_removeClones;
-  bool m_dumpAllHits;
-  bool m_Occupancy;
-  std::string m_hitManagerName;
-  bool m_doClusterResidual;
-  bool m_doTrackStudy;
-  
-  //from PrCounter
-  // typedef LinkerTool<LHCb::Track, LHCb::MCParticle> MyAsct;  
-  // typedef MyAsct::DirectType            Table;  
-  // typedef MyAsct::DirectType::Range     Range;
-  // typedef Table::iterator               iterator;                     
-  // typedef MyAsct::InverseType           InvTable;
-  // typedef InvTable::Range               InvRange;
-  // typedef InvTable::iterator            InvIterator;
-  // MyAsct*         m_link;
-  // const InvTable* m_invTable;                                  
+  //need the ftDet to produce the pseudochannel conversion numbering in occupancy plot for v61
+  DeFTDetector * m_ftDet = nullptr;
 };
 #endif // PRCLUSTERSRESIDUAL_H
diff --git a/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp b/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp
index d7bf33cc1ec4d9e84208f6d7cc8764b3f5c948d3..9b979457bf71e5e91bee30e32dc005cf613ee132 100644
--- a/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp
+++ b/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp
@@ -218,7 +218,7 @@ StatusCode PrDebugTrackingLosses::execute() {
           } else if ( (*itId).isFT() ) {
             LHCb::FTChannelID ftID = (*itId).ftID();
             info() << format( "    FT St%2d La%2d Pm%2d Cel%4d    ",
-                              ftID.layer()/4, ftID.layer()&3, ftID.sipmId(), ftID.sipmCell() );
+                              ftID.station(), ftID.layer(), ftID.sipm(), ftID.channel() );
           } else if ( (*itId).isOT() ) {
             LHCb::OTChannelID otID = (*itId).otID();
             info() << format( "    OT St%2d La%2d mo%2d Str%4d    ",
@@ -285,3 +285,4 @@ StatusCode PrDebugTrackingLosses::finalize() {
 }
 
 //=============================================================================
+
diff --git a/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp b/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp
index ef1d3220cfbd53cca48e54eaec2a99b0552e548b..a7dc95b7e0c07b86a3b9c67f14803f8979b5d8da 100755
--- a/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp
+++ b/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp
@@ -8,7 +8,6 @@
 #include "Event/VPCluster.h"
 #include "Event/STCluster.h"
 #include "Event/OTTime.h"
-#include "Event/FTCluster.h"
 #include "Event/FTLiteCluster.h"
 
 // Range V3
@@ -137,7 +136,7 @@ StatusCode PrLHCbID2MCParticle::execute() {
   //== FT
   const auto* ft_clusters = getIfExists<FastClusterContainer<LHCb::FTLiteCluster,int>>( LHCb::FTLiteClusterLocation::Default );
   if ( ft_clusters ) {
-    LinkedTo<LHCb::MCParticle> ftLink( evtSvc(), msgSvc(),LHCb::FTClusterLocation::Default );
+    LinkedTo<LHCb::MCParticle> ftLink( evtSvc(), msgSvc(),LHCb::FTLiteClusterLocation::Default );
     if ( !ftLink.notFound() ) {
       for(const auto& cluster : *ft_clusters) {
         linkAll( ftLink, lhcbLink, cluster.channelID(), {cluster.channelID()});
diff --git a/Pr/PrMCTools/src/PrPlotFTHits.cpp b/Pr/PrMCTools/src/PrPlotFTHits.cpp
index faff722e190cdc4ae22a4ad3a1ac33a07b90b9f0..974d433549740622a9f8d2e9f9145caa6abf54e2 100644
--- a/Pr/PrMCTools/src/PrPlotFTHits.cpp
+++ b/Pr/PrMCTools/src/PrPlotFTHits.cpp
@@ -114,17 +114,17 @@ void PrPlotFTHits::plotOccupancy(){
     noSpill = true;
   }
   //"/Event/Prev/MC/FT/Hits"
-  LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" );
-
-  LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
-  LinkedTo<LHCb::MCHit, LHCb::FTCluster> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits");
-  //"/Event/Prev/MC/FT/Hits"
-  char layerName[100];
-  char Title[100];
-  std::vector<int> nHits(12,0);
-  MCTrackInfo trackInfo( evtSvc(), msgSvc() );
-
-  LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default );
+  LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" );   
+    
+  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
+  LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits");
+  //"/Event/Prev/MC/FT/Hits"                                             
+  char layerName[100];                                                   
+  char Title[100];                                                       
+  std::vector<int> nHits(12,0);                                          
+  MCTrackInfo trackInfo( evtSvc(), msgSvc() );                           
+    
+  LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default );                                                        
   if(  mcVert == nullptr ){
     error() << "Could not find MCVertices at " << LHCb::MCParticleLocation::Default << endmsg;
   }
@@ -164,19 +164,7 @@ void PrPlotFTHits::plotOccupancy(){
       hit = &*itH;
       nHits[layer]++;
       LHCb::FTLiteCluster liteCluster = getLiteCluster(hit->id());
-      int Module = liteCluster.channelID().module();
-      int Quarter = liteCluster.channelID().quarter();
-      int SipmID = liteCluster.channelID().sipmId();
-      int layer = liteCluster.channelID().layer();
-      int module2 = Module;
-      if(Module == 11 && Module > 4 ) module2 = 1;
-      if(Module != 11 && Module > 4 ) module2 = Module-3;
-      if(Module <  5  ) module2 = 6-Module;
-      if(Module == 10 ) module2 = 1;
-      int index = Quarter*128 + 16*module2+ std::floor(SipmID);
-      // if(index<0 || index>){
-      //   always()<<"Problem in Sipm Index; Correct for that"<<endmsg;
-      // }
+      int index = liteCluster.channelID().uniqueSiPM() & 511;
       sprintf(layerName,"Layer%i/AllContributions",layer);
       sprintf(Title,"Sipm Index Layer %i;SiPm Id;Clusters/Sipm/Event",layer);
       m_histoTool->plot1D(index,layerName,Title,0.5,512.5,512);
@@ -352,7 +340,7 @@ void PrPlotFTHits::plotHitEfficiency(){
   // -- Link different track containers to MCParticles
   LinkedTo<LHCb::MCParticle, LHCb::Track> myForwardLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Forward);
   LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
-  LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
+  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
   MCTrackInfo trackInfo( evtSvc(), msgSvc() );
 
   // -----------------------------------------------------------------------------------
@@ -739,8 +727,8 @@ void PrPlotFTHits::plotMCHits () {
 
 
   // -- MC linking
-  LinkedTo<LHCb::MCHit, LHCb::FTCluster> myMCHitLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits");
-  LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
+  LinkedTo<LHCb::MCHit> myMCHitLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits");
+  LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
 
   MCTrackInfo trackInfo( evtSvc(), msgSvc() );
 
diff --git a/Pr/PrMCTools/src/PrTStationDebugTool.cpp b/Pr/PrMCTools/src/PrTStationDebugTool.cpp
index 304fbea5b1c195731ea92e7fe7f7c34cb29660df..67a93429c5be0e6b59f122e5f4e7a043016765e3 100644
--- a/Pr/PrMCTools/src/PrTStationDebugTool.cpp
+++ b/Pr/PrMCTools/src/PrTStationDebugTool.cpp
@@ -5,7 +5,7 @@
 #include "Linker/LinkedTo.h"
 #include "Kernel/LHCbID.h"
 #include "Event/OTTime.h"
-#include "Event/FTCluster.h"
+#include "Event/FTLiteCluster.h"
 #include "Event/MCParticle.h"
 
 // local
@@ -50,8 +50,8 @@ bool PrTStationDebugTool::matchKey ( LHCb::LHCbID id, int key ) {
       if ( key == part->key() ) return true;
       part = oLink.next();
     }
-  } else if ( id.isFT() ) {
-    LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
+  } else if ( id.isFT() ) {    
+    LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
     LHCb::FTChannelID idO = id.ftID();
 
     LHCb::MCParticle* part = fLink.first( idO );
@@ -77,7 +77,7 @@ void PrTStationDebugTool::printKey( MsgStream& msg, LHCb::LHCbID id ) {
       part = oLink.next();
     }
   } else if ( id.isFT() ) {
-    LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
+    LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
     LHCb::FTChannelID idO = id.ftID();
 
     LHCb::MCParticle* part = fLink.first( idO );
diff --git a/Rec/GlobalReco/root/RichCorrelations.C b/Rec/GlobalReco/root/RichCorrelations.C
new file mode 100644
index 0000000000000000000000000000000000000000..96d0e6757597b3b07360a42106c1b858302f2851
--- /dev/null
+++ b/Rec/GlobalReco/root/RichCorrelations.C
@@ -0,0 +1,66 @@
+
+#include <memory>
+
+void RichCorrelations()
+{
+
+  auto c = std::make_unique<TCanvas>( "RichCor", "RichCor", 1400, 1000 );
+
+  TTree tree;
+  
+  tree.AddFriend("hlt2 = ChargedProtoTuple/protoPtuple",
+                 "/usera/jonesc/LHCbCMake/Brunel/output/future/Quartic100R1R2-AllHypos/protoparticles.tuples.root");
+
+  //tree.AddFriend("hlt1 = ChargedProtoTuple/protoPtuple", "/usera/jonesc/LHCbCMake/Brunel/output/future/FastQuartic50R1R2-AllHypos-NoUnAmBigPhot-NoBeamPipeCheck-1It/protoparticles.tuples.root");
+  //tree.AddFriend("hlt1 = ChargedProtoTuple/protoPtuple", "/usera/jonesc/LHCbCMake/Brunel/output/future/CKEsti50R1R2-KaPi/protoparticles.tuples.root");
+  tree.AddFriend("hlt1 = ChargedProtoTuple/protoPtuple", "/usera/jonesc/LHCbCMake/Brunel/output/future/CKEsti50R1R2-AllHypos/protoparticles.tuples.root");
+
+  tree.Draw("hlt2.TrackP:hlt1.TrackP","hlt2.TrackP<100000","zcol");
+  c->SaveAs( "P-Correlations.pdf" );
+
+  auto sel = TCut("hlt2.RichAbovePiThres>0 && hlt1.RichAbovePiThres>0 && hlt2.TrackP<100000 && hlt2.TrackP>3000");
+
+  auto trueK = TCut("abs(hlt2.MCParticleType) == 321");
+
+  auto truePi = TCut("abs(hlt2.MCParticleType) == 211");
+
+  tree.Draw( "hlt2.RichDLLk:hlt1.RichDLLk",sel && trueK,"zcol");
+  c->SaveAs( "TrueK-RichDLLk-2dCorrelations.pdf" );
+
+  tree.Draw( "hlt2.RichDLLk-hlt1.RichDLLk",sel && trueK );
+  c->SaveAs( "TrueK-RichDLLk-1dCorrelations.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel );
+  c->SaveAs( "All-RichDLLk-HLT1.pdf" );
+
+  tree.Draw( "hlt2.RichDLLk", sel );
+  c->SaveAs( "All-RichDLLk-HLT2.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel && trueK );
+  c->SaveAs( "TrueK-RichDLLk-HLT1.pdf" );
+
+  tree.Draw( "hlt2.RichDLLk", sel && trueK );
+  c->SaveAs( "TrueK-RichDLLk-HLT2.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel && truePi );
+  c->SaveAs( "TruePi-RichDLLk-HLT1.pdf" );
+
+  tree.Draw( "hlt2.RichDLLk", sel && truePi );
+  c->SaveAs( "TruePi-RichDLLk-HLT2.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>-2") );
+  c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt-2.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>0") );
+  c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt0.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>2") );
+  c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt2.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>4") );
+  c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt4.pdf" );
+
+  tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>10") );
+  c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt10.pdf" );
+
+}
diff --git a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C
index 0d62054025282aa02a5d089761c48f3c21f8dbfd..abd6ca3a8a1b9a5e0e1e27f70a5d0edcd8443fcc 100755
--- a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C
+++ b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C
@@ -51,14 +51,14 @@ void RichKaonIDCompareFiles()
 
   //const std::string dir = "/Users/chris/LHCb/RootFiles/Run2";
   //const std::string dir = "/usera/jonesc/LHCbCMake/lhcb-head/BrunelDevNightly/output/Run2";
-  const std::string dir = "/usera/jonesc/LHCbCMake/Brunel/output/future";
+  const std::string dir = "/home/chris/LHCb/future";
 
   typedef std::vector< std::tuple<std::string,std::string,Color_t> > PlotData;
 
   const PlotData plotdata = 
     {
-      std::make_tuple ( "rich", "RICH Current", kBlack  ),
-      std::make_tuple ( "richfuture", "RICH Future", kRed-6  )
+      std::make_tuple ( "old-quartic", "Old Quartic", kBlack  ),
+      std::make_tuple ( "new-quartic", "New NR Quartic", kRed-6  )
     };
 
   // const PlotData plotdata = 
diff --git a/RecSys/CMakeLists.txt b/RecSys/CMakeLists.txt
index 939db4b274ec0abcb564fb5bc7974532298c4d12..fce23855701a083474113c53348326665416e13c 100644
--- a/RecSys/CMakeLists.txt
+++ b/RecSys/CMakeLists.txt
@@ -60,6 +60,7 @@ gaudi_depends_on_subdirs(Calo/CaloMoniDst
                          Rich/RichFutureRecSys
                          Rich/RichFutureRecMonitors
                          Rich/RichFutureRecCheckers
+                         Rich/RichRecTests
                          Tf/FastPV
                          Tf/FastVelo
                          Tf/PatAlgorithms
diff --git a/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py b/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py
index a051222bc35619bb9125deb68eb58812e1aa55a5..af38d34c09b81b0ed0d2d77f80439029c7e3e7ed 100644
--- a/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py
+++ b/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py
@@ -31,7 +31,7 @@ importOptions("$L0TCK/L0DUConfig.opts")
 # For 2016 data
 importOptions("$APPCONFIGOPTS/Brunel/DataType-2016.py")
 LHCbApp().DDDBtag    = "dddb-20150724"
-LHCbApp().CondDBtag  = "cond-20160517"
+LHCbApp().CondDBtag  = "cond-20170120-1"
 
 # Aerogel Sub Tiles
 #CondDB().LocalTags["LHCBCOND"] = ["rich1-20110624"]
diff --git a/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py b/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py
index 516bc597b7f174d8105309b277f845068c01b294..ef223eb1d1e8c9a3b2b8e52863cff120317466aa 100755
--- a/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py
+++ b/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py
@@ -42,6 +42,8 @@ elif 2013 == year :
   ConfigV = ['Collision13','Protonion13','Ionproton13','Ionsmog']
 elif 2015 == year :
   ConfigV = ['Collision15','Collision15em','Protonhelium15','Protonneon15','Protonargon15','Lead15']
+elif 2016 == year :
+  ConfigV = ['Collision16']
 else:
   print 'Unknown year', year
   DIRAC.exit(2)
@@ -51,6 +53,10 @@ allOK = True
 # Dictionary for LFNs by Run
 RunLFNs = { }
 
+polarity = "MagUp"
+#polarity = "MagDown"
+#polarity = "All"
+
 startDate = getDate(year,startmonth,startday)
 endDate   = getDate(year,lastmonth,lastday)
 dict = { 'StartDate'        : startDate,
@@ -81,80 +87,106 @@ else:
   for run in runs :
 
     if maxrun == -1 or run <= maxrun :
-      
+        
       nRun += 1
 
-      for config in ConfigV :
-
-        type = 91000000 # EXPRESS Stream
-        if year == 2009 : type = 90000000 # Use full stream for 2009
-        if run > 71473  and run < 72332 : # Use FULL for first 14nb-1 in 2010
-          type = 90000000 # FULL Stream
-        if run > 77595  and run < 77624 : # Express disappeared for unknown reasons
-          type = 90000000 # FULL Stream
-        if run > 100256 and run < 102177 : # Express turned off by accident after 09/2011 TS
-          type = 90000000 # FULL Stream
-        if config == 'Collision11_25' : # No express for 2011 25ns tests
-          type = 90000000 # FULL Stream
-        if config == 'Protonion12' : # (Currently) no express stream for pA data
-          type = 90000000 # FULL Stream
-        if config == 'Protonion13' or config == 'Ionproton13' :
-          type = 90000000 # FULL Stream
-        #if  run > 77595  and run < 77624 : # Express disappeared for unknown reasons
-        if run > 157596 : # Express stream turned off in Run II
-          type = 90000000 # FULL Stream
-
-        typeS = "EXPRESS"
-        if type == 90000000 : typeS = "FULL"
-
-        # Processing pass. For 2011 express stream data used a Merged location
-        procpass = 'Real Data'
-        if year == 2011 and type != 90000000 :
-          procpass = 'Real Data/Merging'
-
-        # Raw files
-        bkDict = { 'ConfigName'        : 'LHCb',
-                   'ConfigVersion'     : config,
-                   'ProcessingPass'    : procpass,
-                   'FileType'          : 'ALL',
-                   'StartRun'          : run,
-                   'EndRun'            : run,
-                   'EventType'         : type }
-        resultB = database.getFilesWithGivenDataSets(bkDict)
-
-        if not resultB['OK']:
-          print "Error querying BK with", bkDict, resultB['Message']
-          allOK = False
-        else:
-          tmpLFNList = [ ]
-          for lfn in resultB['Value']:
-            filetype = lfn.split('.')[1]
-            if filetype == 'dst' or filetype == 'raw':
-              tmpLFNList += ["LFN:"+lfn.encode("ascii")]
-          tmpLFNList.sort()
-          if len(tmpLFNList) > 0 :
-            if run not in RunLFNs.keys() : RunLFNs[run] = [ ]
-            RunLFNs[run] += tmpLFNList
-
-      if not allOK : break
+      # Polarity ?
+      res = database.getRunInformations(int(run))
+      if polarity == "All" or polarity in res['Value']['DataTakingDescription'] :
+        
+        for config in ConfigV :
+
+          type = 91000000 # EXPRESS Stream
+          if year == 2009 : type = 90000000 # Use full stream for 2009
+          if run > 71473  and run < 72332 : # Use FULL for first 14nb-1 in 2010
+            type = 90000000 # FULL Stream
+          if run > 77595  and run < 77624 : # Express disappeared for unknown reasons
+            type = 90000000 # FULL Stream
+          if run > 100256 and run < 102177 : # Express turned off by accident after 09/2011 TS
+            type = 90000000 # FULL Stream
+          if config == 'Collision11_25' : # No express for 2011 25ns tests
+            type = 90000000 # FULL Stream
+          if config == 'Protonion12' : # (Currently) no express stream for pA data
+            type = 90000000 # FULL Stream
+          if config == 'Protonion13' or config == 'Ionproton13' :
+            type = 90000000 # FULL Stream
+          #if  run > 77595  and run < 77624 : # Express disappeared for unknown reasons
+          if run > 157596 : # Express stream turned off in Run II
+            type = 90000000 # FULL Stream
+          if  config == 'Collision16' :
+            type = 90000000 # FULL Stream
+
+          typeS = "EXPRESS"
+          if type == 90000000 : typeS = "FULL"
+
+          # Processing pass. For 2011 express stream data used a Merged location
+          procpass = 'Real Data'
+          if year == 2011 and type != 90000000 :
+            procpass = 'Real Data/Merging'
+
+          # Raw files
+          bkDict = { 'ConfigName'        : 'LHCb',
+                     'ConfigVersion'     : config,
+                     'ProcessingPass'    : procpass,
+                     'FileType'          : 'ALL',
+                     'StartRun'          : run,
+                     'EndRun'            : run,
+                     'EventType'         : type }
+          resultB = database.getFilesWithGivenDataSets(bkDict)
+
+          if not resultB['OK']:
+
+            print "Error querying BK with", bkDict, resultB['Message']
+            allOK = False
+
+          else:
+
+            tmpLFNList = [ ]
             
-      nLFNs = 0
-      if run in RunLFNs.keys() : nLFNs = len(RunLFNs[run])
-      if nLFNs > 0 :
-        print "(", nRun, "of", len(runs), ") Selected", nLFNs, "LFNs for run", run, ConfigV
-        print RunLFNs[run]
-      else:
-        print "(", nRun, "of", len(runs), ") No data selected for run", run, ConfigV
+            lfns = resultB['Value']
+            tmplfns = list( set(lfns) )
+            if len(tmplfns) != len(lfns) :
+              print "WARNING : run", run, "has duplicate LFNs"
+              print bkDict
+              print resultB['Value']
+              lfns = tmplfns
+
+            for lfn in lfns:
+              filetype = lfn.split('.')[1]
+              if filetype == 'dst' or filetype == 'raw':
+                tmpLFNList += ["LFN:"+lfn.encode("ascii")]
+
+            tmpLFNList.sort()
+            if len(tmpLFNList) > 0 :
+              if run not in RunLFNs.keys() : RunLFNs[run] = [ ]
+              RunLFNs[run] += tmpLFNList
+
+        if not allOK : break
+      
+        nLFNs = 0
+        if run in RunLFNs.keys() : nLFNs = len(RunLFNs[run])
+        if nLFNs > 0 :
+          print "(", nRun, "of", len(runs), ") Selected", nLFNs, "LFNs for run", run, ConfigV
+           #print RunLFNs[run]
+        else:
+          print "(", nRun, "of", len(runs), ") No data selected for run", run, ConfigV
 
-    else:
+      else :
+
+        print "(", nRun, "of", len(runs), ") Run", run, "has wrong polarity"
 
+    else:
+      
       print "Skipping run", run
         
 if allOK :
     
   # Pickle the data
   if len(RunLFNs.keys()) > 0:
-    filename = "RunData/RunLFNs_"+startDate+"_"+endDate+".pck.bz2"
+    if polarity == "All" :
+      filename = "RunData/RunLFNs_"+startDate+"_"+endDate+".pck.bz2"
+    else :
+      filename = "RunData/RunLFNs_"+polarity+"_"+startDate+"_"+endDate+".pck.bz2"
     file = bz2.BZ2File(filename,"w")
     pickle.dump(RunLFNs,file)
     file.close()
diff --git a/Rich/RichFutureGlobalPID/src/RichGlobalPIDNamespaces.h b/Rich/RichFutureGlobalPID/src/RichGlobalPIDNamespaces.h
new file mode 100755
index 0000000000000000000000000000000000000000..d28f230cf03c26afe94d0b2a15302dd4491553c1
--- /dev/null
+++ b/Rich/RichFutureGlobalPID/src/RichGlobalPIDNamespaces.h
@@ -0,0 +1,20 @@
+
+//-----------------------------------------------------------------------------
+/** @namespace Rich::Future::Rec::GlobalPID
+ *
+ *  Namespace for RICH Global PID software
+ *
+ *  @author Chris Jones  Christopher.Rob.Jones@cern.ch
+ *  @date   04/12/2006
+ */
+//-----------------------------------------------------------------------------
+
+//-----------------------------------------------------------------------------
+/** @namespace Rich::Future::Rec::GlobalPID::MC
+ *
+ *  Namespace for RICH Global PID MC related software
+ *
+ *  @author Chris Jones  Christopher.Rob.Jones@cern.ch
+ *  @date   05/12/2006
+ */
+//-----------------------------------------------------------------------------
diff --git a/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h b/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h
index c7472cd147a6058b29182d4282d99fc7b9d91dc8..8443cac3c582fd1d01833b22cc7ce3059a0fd789 100644
--- a/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h
+++ b/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h
@@ -5,7 +5,7 @@
  *  General namespace for all RICH reconstruction software
  *
  *  @author Chris Jones  Christopher.Rob.Jones@cern.ch
- *  @date   08/07/2004
+ *  @date   01/02/2017
  */
 //-----------------------------------------------------------------------------
 
@@ -15,6 +15,16 @@
  *  General namespace for RICH reconstruction MC related software
  *
  *  @author Chris Jones  Christopher.Rob.Jones@cern.ch
- *  @date   05/12/2006
+ *  @date   01/02/2017
+ */
+//-----------------------------------------------------------------------------
+
+//-----------------------------------------------------------------------------
+/** @namespace Rich::Future::Rec::Moni
+ *
+ *  General namespace for RICH reconstruction monitoring utilities
+ *
+ *  @author Chris Jones  Christopher.Rob.Jones@cern.ch
+ *  @date   01/02/2017
  */
 //-----------------------------------------------------------------------------
diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h
index 45351b5eed6973d666964a669290f0abf8d9f578..d7391683b51c73dc3398edb6b5c3b9f741c4fc8c 100644
--- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h
+++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h
@@ -24,6 +24,7 @@ namespace Rich
       /// Type for flags to say if a photon is active or not. i.e. has any signal.
       using PhotonFlags = LHCb::STL::Vector<bool>;
 
+      /// photon active flags TES locations
       namespace PhotonActiveFlagsLocation
       {
         /// Location in TES for the active photon flags
diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h
index 0c5854a24002c19c9d34e5bf6355306f563f6ce9..2c322afa970636a85e300d04a42c2efb5f72200b 100644
--- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h
+++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h
@@ -24,6 +24,14 @@ namespace Rich
     { 
       // ===========================================================================
       
+      /** @namespace Summary
+       *
+       *  General namespace for RICH reconstruction summary information
+       *
+       *  @author Chris Jones  Christopher.Rob.Jones@cern.ch
+       *  @date   01/02/2017
+       */
+
       namespace Summary
       {
         
@@ -125,6 +133,14 @@ namespace Rich
         
         // ===========================================================================
         
+        /** @class Pixel RichSummaryEventData.h
+         *
+         *  Summary of the immutable pixell data.
+         *
+         *  @author Chris Jones
+         *  @date   2016-09-30
+         */
+
         class Pixel final
         {
           
@@ -155,7 +171,7 @@ namespace Rich
         
         // ===========================================================================
         
-        /// TES locations
+        /// Summary information TES locations
         namespace TESLocations
         {
           /// Tracks
diff --git a/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt b/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt
index 5ccf26d32ef7c847a349b75e505bf3d92e17bd4d..85678715d1a9163c5a86e75d796e0601aea34a63 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt
+++ b/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt
@@ -5,7 +5,8 @@ gaudi_subdir(RichFutureRecPhotonAlgorithms v1r0)
 
 gaudi_depends_on_subdirs(Rich/RichFutureRecBase
                          Rich/RichFutureRecEvent
-                         Rich/RichUtils)
+                         Rich/RichUtils
+                         Kernel/VectorClass)
 
 find_package(Boost)
 find_package(GSL)
@@ -16,10 +17,23 @@ find_package(ROOT)
 include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS}
                                                  ${EIGEN_INCLUDE_DIRS})
 
-# For testing high vectorisation and optimisation levels
-#set_property(SOURCE src/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -msse4 -O3 " )
-
 gaudi_add_module(RichFutureRecPhotonAlgorithms
-                 src/*.cpp
-                 INCLUDE_DIRS Boost GSL VDT Eigen Rich/RichFutureRecBase Rich/RichFutureRecEvent Rich/RichUtils
+                 src/*.cpp src/generic/*.cpp src/sse4/*.cpp src/avx/*.cpp src/avx2/*.cpp
+                 INCLUDE_DIRS Boost GSL VDT Eigen Kernel/VectorClass Rich/RichFutureRecBase Rich/RichFutureRecEvent Rich/RichUtils
                  LINK_LIBRARIES Boost GSL VDT RichFutureRecBase RichFutureRecEvent RichUtils)
+
+set_property(SOURCE src/generic/RichQuarticPhotonReco.cpp  APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes " )
+
+set_property(SOURCE src/sse4/RichQuarticPhotonReco.cpp     APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -msse4.2 " )
+
+set_property(SOURCE src/avx/RichQuarticPhotonReco.cpp      APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -mavx " )
+
+# Until such a time that all build targets support AVX2+FMA, fallback to just AVX here
+#set_property(SOURCE src/avx2/RichQuarticPhotonReco.cpp     APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -mavx2 -mfma " )
+set_property(SOURCE src/avx2/RichQuarticPhotonReco.cpp     APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -mavx " )
+
+if(LCG_COMP STREQUAL "gcc" OR
+    (BINARY_TAG_COMP_NAME STREQUAL "gcc" AND BINARY_TAG_COMP_VERSION VERSION_LESS "5"))
+  set_property(SOURCE src/avx/RichQuarticPhotonReco.cpp  APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " )
+  set_property(SOURCE src/avx2/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " )
+endif()
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h
index ff0522a36413bd0cf966d4cc76528d816e0c6b94..09fadbbfe24f72ed61239dc1587e190d6cc7aeeb 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h
@@ -166,18 +166,18 @@ namespace Rich
         }
         
         /// Check the final Cherenkov angles
-        template < typename HYPODATA >
+        template < typename HYPODATA, typename FTYPE >
         inline bool checkAngles( const Rich::RadiatorType rad,
                                  const HYPODATA& tkCkAngles,
                                  const HYPODATA& tkCkRes,
-                                 const float ckTheta,
-                                 const float ckPhi ) const noexcept
+                                 const FTYPE ckTheta,
+                                 const FTYPE ckPhi ) const noexcept
         {
           return (
             // First the basic checks
-            ckPhi   > 0 &&
-            ckTheta > absMinCKTheta(rad) &&
             ckTheta < absMaxCKTheta(rad) &&
+            ckTheta > absMinCKTheta(rad) &&
+            ckPhi   > 0                  &&
             // Now check each hypo
             std::any_of( activeParticles().begin(),
                          activeParticles().end(),
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
index 595911486692d9fc8fdac5ed3215862d017de142..14f6b40195c42be4b526f039f9ff8c0c58e7dc09 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
@@ -47,6 +47,18 @@ StatusCode QuarticPhotonReco::initialize()
   auto sc = MultiTransformer::initialize();
   if ( !sc ) return sc;
 
+  // Get the CPU capabilities and set dispatch method
+  if ( m_detectCPU )
+  {
+    enum CPUID { GENERIC = 0, SSE4 = 6, AVX = 7, AVX2 = 8 };
+    const auto cpuLevel = instrset_detect();
+    _ri_debug << "Instruction set level = " << cpuLevel << endmsg;
+    if      ( cpuLevel >= AVX2 ) { m_run = &QuarticPhotonReco::run_avx2; }
+    else if ( cpuLevel >= AVX  ) { m_run = &QuarticPhotonReco::run_avx; }
+    else if ( cpuLevel >= SSE4 ) { m_run = &QuarticPhotonReco::run_sse4; }
+    else                         { m_run = &QuarticPhotonReco::run_generic; }
+  }
+
   // get the detector elements
   m_rich[Rich::Rich1] = getDet<DeRich>( DeRichLocations::Rich1 );
   m_rich[Rich::Rich2] = getDet<DeRich>( DeRichLocations::Rich2 );
@@ -112,678 +124,6 @@ StatusCode QuarticPhotonReco::initialize()
 
 //=============================================================================
 
-OutData
-QuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segments,
-                               const CherenkovAngles::Vector& ckAngles,
-                               const CherenkovResolutions::Vector& ckResolutions,
-                               const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
-                               const SegmentPhotonFlags::Vector& segPhotFlags,
-                               const Rich::PDPixelCluster::Vector& clusters,
-                               const SpacePointVector& globalPoints,
-                               const SpacePointVector& localPoints,
-                               const Relations::TrackToSegments::Vector& tkToSegRels ) const
-{
-  // use ranges v3 library
-  using namespace ranges::v3;
-
-  // make the output data
-  OutData outData;
-
-  // Shortcut to the photons and relations
-  auto & photons   = std::get<0>(outData);
-  auto & relations = std::get<1>(outData);
-  // guess at reserve size
-  const auto resSize = segments.size() * globalPoints.size() / 20;
-  photons.reserve( resSize );
-  relations.reserve( resSize );
-
-  // View range for pixel data
-  const auto allPixels = Ranges::ConstZip( clusters, globalPoints, localPoints );
-
-  // Make 'ranges' for RICH pixels
-  const auto RPixs = richRanges( allPixels );
-  const auto & RichRanges  = RPixs[0];   // Complete RICH ranges
-  const auto & Rich1Ranges = RPixs[1];   // ranges for RICH1 [top,bottom]
-  const auto & Rich2Ranges = RPixs[2];   // ranges for RICH2 [left,right]
-  
-  // local position corrector
-  // longer term need to remove this
-  const Rich::Rec::RadPositionCorrector corrector;
-
-  // global photon index
-  int photonIndex(-1);
-  
-  // Loop over the track->segment relations
-  for ( const auto & inTkRel : tkToSegRels )
-  {
-    // loop over segments for this track
-    for ( const auto & segIndex : inTkRel.segmentIndices )
-    {
-      // Get the data from the segment containers
-      const auto & segment     = segments[segIndex];
-      const auto & tkCkAngles  = ckAngles[segIndex];
-      const auto & tkCkRes     = ckResolutions[segIndex];
-      const auto & tkLocPtn    = trHitPntsLoc[segIndex];
-      const auto & segFlags    = segPhotFlags[segIndex];
-
-      // Is this segment above threshold. Should be implicit due to above segment loop
-      //if ( ! ( tkCkAngles[lightestActiveHypo()] > 0 ) ) continue;
-
-      // which RICH and radiator
-      const auto rich = segment.rich();
-      const auto rad  = segment.radiator();
-
-      // This implementiation does not support Aerogel
-      if ( UNLIKELY( Rich::Aerogel == rad ) ) { Exception("Aerogel not supported"); }
-
-      // get the best pixel range for this segment, based on where hits are expected
-      const auto& pixR = ( segFlags.inBothPanels() ? RichRanges[rich] :
-                           Rich::Rich1 == rich ? 
-                           ( segFlags.inPanel(Rich::top) ? 
-                             Rich1Ranges[Rich::top]  : Rich1Ranges[Rich::bottom] ) :
-                           ( segFlags.inPanel(Rich::left) ? 
-                             Rich2Ranges[Rich::left] : Rich2Ranges[Rich::right] ) );
-
-      // Loop over pixel information
-      int pixIndex = pixR.first - allPixels.begin() - 1; // (start index - 1) for this range
-      for ( auto pix_it = pixR.first; pix_it != pixR.second; ++pix_it )
-      {
-        // load the data from the tuple
-        const auto & cluster = std::get<0>(*pix_it);
-        const auto & gloPos  = std::get<1>(*pix_it);
-        const auto & locPos  = std::get<2>(*pix_it);
-        // increment the pixel index
-        ++pixIndex;
-
-        //_ri_verbo << "Trying Segment P=" << segment.bestMomentum() << " Pixel " 
-        //          << cluster << endmsg;
-
-        // Apply pre-selection
-      
-        // Pixel position, in local HPD coords corrected for average radiator distortion
-        // Might need to pre-compute and cache this as we repeat it here each segment...
-        const auto pixP = corrector.correct(locPos,rad);
-
-        //  Track local hit point on the same panel as the hit
-        /// @todo CRJ Need to figure out why I was not correcting this in the current stack...
-        ///           Probably correct but should check. 
-        //const auto segPanelPnt = corrector.correct(tkLocPtn.point(cluster.panel()),rad);
-        const auto& segPanelPnt = tkLocPtn.point(cluster.panel());
-
-        // compute the seperation squared
-        const auto sep2 = ( std::pow( (pixP.x()-segPanelPnt.x()), 2 ) +
-                            std::pow( (pixP.y()-segPanelPnt.y()), 2 ) );
-
-        // presel. default to rejected
-        bool SelOK = false;
-
-        // Check overall boundaries
-        //_ri_verbo << " -> sep2 = " << sep2 << " " << m_maxROI2PreSel[rad]
-        //          << " " << m_minROI2PreSel[rad] << endmsg;
-        if ( sep2 < m_maxROI2PreSel[rad] && sep2 > m_minROI2PreSel[rad] )
-        {
-        
-          // estimated CK theta
-          const auto ckThetaEsti = std::sqrt(sep2)*m_scalePreSel[rad];
-          //_ri_verbo << " -> CK theta Esti " << ckThetaEsti << endmsg;
-
-          // Is the hit close to any mass hypo in local coordinate space
-          SelOK = any_of( activeParticles(),
-                         [sigma=m_nSigmaPreSel[rad],&ckThetaEsti,&tkCkAngles,&tkCkRes]
-                         ( const auto hypo )
-                         { return fabs(tkCkAngles[hypo]-ckThetaEsti) < sigma*tkCkRes[hypo]; } );
-          
-        } // boundary check
-
-        // Did we pass the pre-sel
-        if ( !SelOK )
-        {
-          //_ri_verbo << "   -> FAILED pre-selection -> reject" << endmsg;
-          continue;
-        }
-
-        // Move on to the full Quartic reconstruction
-
-        //std::cout << "Starting photon reco" << std::endl;
-
-        // make a photon object to work on
-        photons.emplace_back();
-        auto & gPhoton = photons.back();
-        DataGuard<decltype(photons)> photGuard(photons);
-
-        // pixel side
-        const auto side = cluster.panel();
-
-        // Emission point to use for photon reconstruction
-        // operate directly on photon data member for efficiency
-        //auto & emissionPoint = gPhoton.emissionPoint();
-        //emissionPoint = segment.bestPoint();
-        auto emissionPoint = segment.bestPoint();
-
-        // Photon direction at emission point
-        // Again, operator directly on data member
-        //auto & photonDirection = gPhoton.emissionDir();
-        Gaudi::XYZVector photonDirection;
-
-        // Final reflection points on sec and spherical mirrors
-        // operate directly on photon data
-        //auto & sphReflPoint = gPhoton.sphMirReflectionPoint();
-        //auto & secReflPoint = gPhoton.flatMirReflectionPoint();
-        Gaudi::XYZPoint sphReflPoint, secReflPoint;
-
-        // fraction of segment path length accessible to the photon
-        float fraction(1);
-
-        // Pointers to best sec and spherical mirror segments
-        const DeRichSphMirror * sphSegment = nullptr;
-        const DeRichSphMirror * secSegment = nullptr;
-
-        // flag to say if this photon candidate is un-ambiguous - default to false
-        bool unambigPhoton = false;
-
-        // find the reflection of the detection point in the sec mirror
-        // (virtual detection point) starting with nominal values
-        // At this we are assuming a flat nominal mirror common to all segments
-        auto distance = m_rich[rich]->nominalPlane(side).Distance(gloPos);
-        auto virtDetPoint = gloPos - 2.0 * distance * m_rich[rich]->nominalNormal(side);
-
-        // --------------------------------------------------------------------------------------
-        // For gas radiators, try start and end points to see if photon is unambiguous
-        // NOTE : For this we are using the virtual detection point determined using
-        // the noimnal flat secondary mirror plane. Now the secondary mirrors are actually
-        // spherical this may be introducing some additional uncertainties.
-        // --------------------------------------------------------------------------------------
-        if ( m_testForUnambigPhots[rad] )
-        {
-          if ( UNLIKELY( rad == Rich::Aerogel ) )
-          {
-            // use default emission point and assume unambiguous since path length is so short..
-            unambigPhoton = true;
-          }
-          else  // gas radiators
-          {
-          
-            // flag for beam pipe check. Default is OK
-            bool beamTestOK(true);
-          
-            // -------------------------------------------------------------------------------
-            // First emission point, at start of track segment
-            const auto emissionPoint1 = segment.bestPoint(0.01);
-            // Find mirror segments for this emission point
-            const DeRichSphMirror* sphSegment1 = nullptr;
-            const DeRichSphMirror* secSegment1 = nullptr;
-            Gaudi::XYZPoint sphReflPoint1, secReflPoint1;
-            if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint1,
-                                  sphSegment1, secSegment1, sphReflPoint1, secReflPoint1 ) )
-            {
-              //_ri_debug << rad << " : Failed to reconstruct photon for start of segment" << endmsg;
-              continue;
-            }
-            if ( m_checkBeamPipe[rad] )
-            {
-              beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint1, sphReflPoint1 );
-            }
-            // -------------------------------------------------------------------------------
-          
-            // -------------------------------------------------------------------------------
-            // now do it again for emission point #2, at end of segment
-            const auto emissionPoint2 = segment.bestPoint(0.99);
-            // Find mirror segments for this emission point
-            const DeRichSphMirror* sphSegment2 = nullptr;
-            const DeRichSphMirror* secSegment2 = nullptr;
-            Gaudi::XYZPoint sphReflPoint2, secReflPoint2;
-            if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint2,
-                                  sphSegment2, secSegment2, sphReflPoint2, secReflPoint2 ) )
-            {
-              //_ri_debug << rad << " : Failed to reconstruct photon for end of segment" << endmsg;
-              continue;
-            }
-            if ( !beamTestOK && m_checkBeamPipe[rad] )
-            {
-              beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint2, sphReflPoint2 );
-            }
-            // -------------------------------------------------------------------------------
-          
-            // -------------------------------------------------------------------------------
-            if ( !beamTestOK )
-            {
-              // both start and end points failed beam pipe test -> reject
-              //_ri_debug << rad << " : Failed ambiguous photon beampipe intersection checks" << endmsg;
-              continue;
-            }
-            // -------------------------------------------------------------------------------
-          
-            // -------------------------------------------------------------------------------
-            // Test to see if photon direction hits the real physical mirror segments
-            // -------------------------------------------------------------------------------
-            if ( UNLIKELY(m_checkPrimMirrSegs[rad]) )
-            {
-              // primary mirrors
-              const auto ok = ( sphSegment1->intersects( emissionPoint1,
-                                                         sphReflPoint1-emissionPoint1 ) ||
-                                sphSegment2->intersects( emissionPoint2,
-                                                         sphReflPoint2-emissionPoint2 ) );
-              if ( !ok )
-              {
-                //_ri_debug << rad << " : Failed mirror segment intersection checks" << endmsg;
-                continue;
-              }
-            }
-            // -------------------------------------------------------------------------------
-          
-            // -------------------------------------------------------------------------------
-            // Get the best gas emission point
-            if ( !getBestGasEmissionPoint( rad,
-                                           sphReflPoint1,
-                                           sphReflPoint2,
-                                           gloPos,
-                                           segment,
-                                           emissionPoint,
-                                           fraction ) )
-            {
-              //_ri_debug << rad << " : Failed to compute best gas emission point" << endmsg;
-              continue;
-            }
-            // -------------------------------------------------------------------------------
-          
-            // -------------------------------------------------------------------------------
-            // Is this an unambiguous photon - I.e. only one possible mirror combination
-            if ( ( sphSegment1 == sphSegment2 ) && ( secSegment1 == secSegment2 ) )
-            {
-              // Set pointers to the mirror detector objects
-              sphSegment = sphSegment1;
-              secSegment = secSegment1;
-              // rough guesses at reflection points (improved later on)
-              sphReflPoint = sphReflPoint1 + 0.5*(sphReflPoint2-sphReflPoint1);
-              secReflPoint = secReflPoint1 + 0.5*(secReflPoint2-secReflPoint1);
-              // photon is not unambiguous
-              unambigPhoton = true;
-            }
-            // -------------------------------------------------------------------------------
-          
-          } // end radiator type if
-        
-          // if configured to do so reject ambiguous photons
-          if ( UNLIKELY( m_rejectAmbigPhots[rad] && !unambigPhoton ) )
-          {
-            //_ri_debug << rad << " : Failed ambiguous photon test" << endmsg;
-            continue;
-          }
-        
-        } // end unambiguous photon check
-        // --------------------------------------------------------------------------------------
-
-        // --------------------------------------------------------------------------------------
-        // Active segment fraction cut
-        // --------------------------------------------------------------------------------------
-        if ( UNLIKELY( fraction < m_minActiveFrac[rad] ) )
-        {
-          //_ri_debug << rad << " : Failed active segment fraction cut" << endmsg;
-          continue;
-        }
-        // --------------------------------------------------------------------------------------
-
-        // --------------------------------------------------------------------------------------
-        // if ambiguous gas photon or if ambiguous photon check above has been skipped, try again
-        // using best emission point and nominal mirror geometries to get the spherical and sec
-        // mirrors. Also, force this reconstruction if the above unambiguous test was skipped
-        // --------------------------------------------------------------------------------------
-        if ( UNLIKELY( !m_testForUnambigPhots[rad] || !unambigPhoton ) )
-        {
-          if ( !findMirrorData(rich,side,
-                               virtDetPoint,emissionPoint,
-                               sphSegment,secSegment,sphReflPoint,secReflPoint ) )
-          {
-            //_ri_debug << rad << " : Failed backup photon reconstruction" << endmsg;
-            continue;
-          }
-        }
-        // --------------------------------------------------------------------------------------
-      
-        // --------------------------------------------------------------------------------------
-        // Finally reconstruct the photon using best emission point and the best mirror segments
-        // --------------------------------------------------------------------------------------
-        if ( m_useAlignedMirrSegs[rad] )
-        {
-        
-          // If iterations are disabled for this radiator, update the virtual detection point
-          // using the selected mirror segment
-          if ( UNLIKELY( 0 == m_nMaxQits[rad] ) )
-          {
-          
-            // Form the virtual detection point
-            distance     = secSegment->centreNormalPlane().Distance(gloPos);
-            virtDetPoint = gloPos - 2.0 * distance * secSegment->centreNormal();
-          
-            // solve the quartic using the new data
-            m_quarticSolver.solve<double>( emissionPoint,
-                                           sphSegment->centreOfCurvature(),
-                                           virtDetPoint,
-                                           sphSegment->radius(),
-                                           sphReflPoint );
-          
-            // (re)find the spherical mirror segment
-            // CRJ - Is this needed ?
-            sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
-          
-          }
-          else
-          {
-          
-            // Iterate to final solution, improving the secondary mirror info
-            int iIt(0);
-            Gaudi::XYZPoint last_mirror_point(0,0,0);
-            bool loopOK = true;
-            while ( iIt < m_nMaxQits[rad] )
-            {
-            
-              // Get secondary mirror reflection point,
-              // using the best actual secondary mirror segment at this point
-              const auto dir = virtDetPoint - sphReflPoint;
-              const auto sc = ( m_treatSecMirrsFlat[rich] ? 
-                                intersectPlane( sphReflPoint, dir,
-                                                secSegment->centreNormalPlane(),
-                                                secReflPoint ) :
-                                intersectSpherical( sphReflPoint, dir,
-                                                    secSegment->centreOfCurvature(),
-                                                    secSegment->radius(),
-                                                    secReflPoint ) );
-              if ( UNLIKELY(!sc) )
-              {
-                //_ri_debug << rad << " : Failed to intersect nominal secondary mirror plane" << endmsg;
-                loopOK = false;
-                break;
-              }
-
-              // (re)find the secondary mirror
-              secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint );
-            
-              // Construct plane tangential to secondary mirror passing through reflection point
-              const Gaudi::Plane3D plane( secSegment->centreOfCurvature()-secReflPoint, secReflPoint );
-            
-              // re-find the reflection of the detection point in the sec mirror
-              // (virtual detection point) with this mirror plane
-              distance     = plane.Distance(gloPos);
-              virtDetPoint = gloPos - 2.0 * distance * plane.Normal();
-            
-              // solve the quartic using the new data
-              m_quarticSolver.solve<double>( emissionPoint,
-                                             sphSegment->centreOfCurvature(),
-                                             virtDetPoint,
-                                             sphSegment->radius(),
-                                             sphReflPoint );
-
-              // (re)find the spherical mirror segment
-              sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
-
-              // for iterations after the first, see if we are still moving
-              // If not, abort iterations early
-              if ( iIt > m_nMinQits[rad] &&
-                   (last_mirror_point-secReflPoint).Mag2() < m_minSphMirrTolIt[rad] ) { break; }
-
-              // store last sec mirror point
-              last_mirror_point = secReflPoint;
-
-              // increment iteration counter
-              ++iIt;
-            }
-            // if above while loop failed, abort this photon
-            if ( UNLIKELY(!loopOK) ) { continue; }
-          
-          }
-        
-        }
-        // --------------------------------------------------------------------------------------
-      
-        // --------------------------------------------------------------------------------------
-        // check that spherical mirror reflection point is on the same side as detection point
-        // and (if configured to do so) photon does not cross between detector sides
-        // --------------------------------------------------------------------------------------
-        if ( UNLIKELY( !sameSide(rad,sphReflPoint,virtDetPoint) ) )
-        {
-          //_ri_debug << rad << " : Reflection point on wrong side" << endmsg;
-          continue;
-        }
-        if ( UNLIKELY( m_checkPhotCrossSides[rad] &&
-                       !sameSide(rad,sphReflPoint,emissionPoint) ) )
-        {
-          //_ri_debug << rad << " : Photon crosses between detector sides" << endmsg;
-          continue;
-        }
-        //else { continue; } // uncomment to select ONLY crossing photons
-        // --------------------------------------------------------------------------------------
-      
-        // --------------------------------------------------------------------------------------
-        // For as radiators if ambiguous photon checks are disabled (since this is
-        // already done for these photons during those checks), check if the photon would have
-        // intersected with the beampipe
-        // --------------------------------------------------------------------------------------
-        if ( UNLIKELY( m_checkBeamPipe[rad] && !m_testForUnambigPhots[rad] ) )
-        {
-          if ( UNLIKELY( deBeam(rich)->testForIntersection( emissionPoint, sphReflPoint ) ) )
-          {
-            //_ri_debug << rad << " : Failed final beampipe intersection checks" << endmsg;
-            continue;
-          }
-        }
-        // --------------------------------------------------------------------------------------
-
-        // --------------------------------------------------------------------------------------
-        // If using aligned mirror segments, get the final sec mirror reflection
-        // point using the best mirror segments available at this point
-        // For RICH2, use the spherical nature of the scondary mirrors
-        // For RICH1, where they are much flatter, assume complete flatness
-        // --------------------------------------------------------------------------------------
-        if ( m_useAlignedMirrSegs[rad] )
-        {
-          const auto dir = virtDetPoint - sphReflPoint;
-          const auto sc = ( m_treatSecMirrsFlat[rich] ? 
-                            intersectPlane( sphReflPoint, dir,
-                                            secSegment->centreNormalPlane(),
-                                            secReflPoint ) :
-                            intersectSpherical( sphReflPoint, dir,
-                                                secSegment->centreOfCurvature(),
-                                                secSegment->radius(),
-                                                secReflPoint ) );
-          if ( !sc )
-          {
-            //_ri_debug << rad << " : Failed final secondary mirror plane intersection" << endmsg;
-            continue;
-          }
-        }
-        // --------------------------------------------------------------------------------------
-
-        // --------------------------------------------------------------------------------------
-        // calculate the cherenkov angles using the photon and track vectors
-        // --------------------------------------------------------------------------------------
-        photonDirection = (sphReflPoint-emissionPoint).Unit();
-        float thetaCerenkov(0), phiCerenkov(0);
-        segment.angleToDirection( photonDirection, thetaCerenkov, phiCerenkov );
-        // --------------------------------------------------------------------------------------
-
-        //---------------------------------------------------------------------------------------
-        // Apply fudge factor correction for small biases in CK theta
-        //---------------------------------------------------------------------------------------
-        thetaCerenkov += ckThetaCorrection(rad);
-        //---------------------------------------------------------------------------------------
-
-        // --------------------------------------------------------------------------------------
-        // Final checks on the Cherenkov angles
-        // --------------------------------------------------------------------------------------
-        if ( UNLIKELY( !checkAngles( rad, tkCkAngles, tkCkRes,
-                                     thetaCerenkov, phiCerenkov ) ) )
-        {
-          //_ri_verbo << "    -> photon FAILED checkAngleInRange test" << endmsg;
-          continue;
-        }
-        // --------------------------------------------------------------------------------------
-
-        // --------------------------------------------------------------------------------------
-        // Set (remaining) photon parameters
-        // --------------------------------------------------------------------------------------
-        gPhoton.setCherenkovTheta         ( thetaCerenkov       );
-        gPhoton.setCherenkovPhi           ( phiCerenkov         );
-        gPhoton.setActiveSegmentFraction  ( fraction            );
-        //gPhoton.setDetectionPoint         ( gloPos              );
-        gPhoton.setSmartID                ( cluster.primaryID() );
-        gPhoton.setUnambiguousPhoton      ( unambigPhoton       );
-        //gPhoton.setPrimaryMirror          ( sphSegment          );
-        //gPhoton.setSecondaryMirror        ( secSegment          );
-        //_ri_verbo << "Created photon " << gPhoton << endmsg;
-        // --------------------------------------------------------------------------------------
-     
-        // If we get here, keep the just made photon
-        photGuard.setOK();
-        // Save relations
-        relations.emplace_back( ++photonIndex, pixIndex, segIndex, inTkRel.tkIndex );
-        //info() << relations.back() << endmsg;
-
-      }
-
-    }
-    
-  }
-  
-  //_ri_debug << "Created " << photons.size() << " photons" << endmsg;
-  return outData;
-}
-
-//=========================================================================
-// Find mirror segments and reflection points for given data
-//=========================================================================
-bool
-QuarticPhotonReco::
-findMirrorData( const Rich::DetectorType rich,
-                const Rich::Side side,
-                const Gaudi::XYZPoint& virtDetPoint,
-                const Gaudi::XYZPoint& emissionPoint,
-                const DeRichSphMirror*& sphSegment,
-                const DeRichSphMirror*& secSegment,
-                Gaudi::XYZPoint& sphReflPoint,
-                Gaudi::XYZPoint& secReflPoint ) const
-{
-  // solve quartic equation with nominal values and find spherical mirror reflection point
-  m_quarticSolver.solve<double>( emissionPoint,
-                                 m_rich[rich]->nominalCentreOfCurvature(side),
-                                 virtDetPoint,
-                                 m_rich[rich]->sphMirrorRadius(),
-                                 sphReflPoint );
-  // find the spherical mirror segment
-  sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
-  // Search for the secondary segment
-  // Direction vector from primary mirror point to virtual detection point
-  const auto dir ( virtDetPoint - sphReflPoint );
-  // find the sec mirror intersction point and secondary mirror segment
-  const bool OK = intersectPlane( sphReflPoint, dir,
-                                  m_rich[rich]->nominalPlane(side),
-                                  secReflPoint );
-  if ( OK )
-  {
-    // find the secondary mirror
-    secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint );
-    // Re-find the secondary mirror reflection point using the new mirror info
-    // OK = ( !m_treatSecMirrsFlat[rich] ? 
-    //        intersectSpherical( sphReflPoint, dir,
-    //                            secSegment->centreOfCurvature(),
-    //                            secSegment->radius(),
-    //                            secReflPoint ) :
-    //        intersectPlane( sphReflPoint, dir,
-    //                        secSegment->centreNormalPlane(),
-    //                        secReflPoint ) );
-  }
-  // return
-  return OK;
-}
-
-//=========================================================================
-// Compute the best emission point for the gas radiators using
-// the given spherical mirror reflection points
-//=========================================================================
-bool
-QuarticPhotonReco::
-getBestGasEmissionPoint( const Rich::RadiatorType radiator,
-                         const Gaudi::XYZPoint& sphReflPoint1,
-                         const Gaudi::XYZPoint& sphReflPoint2,
-                         const Gaudi::XYZPoint& detectionPoint,
-                         const LHCb::RichTrackSegment & segment,
-                         Gaudi::XYZPoint & emissionPoint,
-                         float & fraction ) const
-{
-  double alongTkFrac = 0.5;
-
-  if ( radiator == Rich::Rich1Gas )
-  {
-    // First reflection and hit point on same y side ?
-    const bool sameSide1 = ( sphReflPoint1.y() * detectionPoint.y() > 0 );
-    const bool sameSide2 = ( sphReflPoint2.y() * detectionPoint.y() > 0 );
-    if ( sameSide1 && sameSide2 )
-    {
-      emissionPoint = segment.bestPoint();
-    }
-    else if ( sameSide1 )
-    {
-      fraction = (float)(std::fabs(sphReflPoint1.y()/(sphReflPoint1.y()-sphReflPoint2.y())));
-      alongTkFrac = fraction/2.0;
-      emissionPoint = segment.bestPoint(alongTkFrac);
-    }
-    else if ( sameSide2 )
-    {
-      fraction = (float)(std::fabs(sphReflPoint2.y()/(sphReflPoint1.y()-sphReflPoint2.y())));
-      alongTkFrac = 1.0-fraction/2.0;
-      emissionPoint = segment.bestPoint(alongTkFrac);
-    }
-    else
-    {
-      //Warning( "Rich1Gas : Both RICH spherical mirror hits opposite side to detection point" );
-      return false;
-    }
-
-  }
-  else if ( radiator == Rich::Rich2Gas )
-  {
-    // First sphReflPoint and hit point on same x side ?
-    const bool sameSide1 = ( sphReflPoint1.x() * detectionPoint.x() > 0 );
-    const bool sameSide2 = ( sphReflPoint2.x() * detectionPoint.x() > 0 );
-    if ( sameSide1 && sameSide2 )
-    {
-      emissionPoint = segment.bestPoint();
-    }
-    else if ( sameSide1 )
-    {
-      fraction = (float)(std::fabs(sphReflPoint1.x()/(sphReflPoint1.x()-sphReflPoint2.x())));
-      alongTkFrac = fraction/2.0;
-      emissionPoint = segment.bestPoint(alongTkFrac);
-    }
-    else if ( sameSide2 )
-    {
-      fraction = (float)(std::fabs(sphReflPoint2.x()/(sphReflPoint1.x()-sphReflPoint2.x())));
-      alongTkFrac = 1.0-fraction/2.0;
-      emissionPoint = segment.bestPoint(alongTkFrac);
-    }
-    else
-    {
-      //Warning( "Rich2Gas : Both RICH spherical mirror hits opposite side to detection point" );
-      return false;
-    }
-  }
-  else { Error( "::getBestGasEmissionPoint() called for Aerogel segment !!" ); }
-
-  // if ( msgLevel(MSG::VERBOSE) )
-  // {
-  //   verbose() << radiator << " best emission point correction :- " << endmsg
-  //             << " -> Photon detection point = " << detectionPoint << endmsg
-  //             << " -> Sph. Mirror ptns       = " << sphReflPoint1 << " " << sphReflPoint2 << endmsg
-  //             << " -> Segment entry/exit     = " << trSeg.entryPoint() << " " << trSeg.exitPoint() << endmsg
-  //             << " -> Segment fraction       = " << fraction << endmsg
-  //             << " -> Emm. Ptn. Along traj   = " << alongTkFrac << endmsg
-  //             << " -> Best Emission point    = " << emissionPoint << endmsg;
-  // }
-
-  return true;
-}
-
-//=============================================================================
-
 // Declaration of the Algorithm Factory
 DECLARE_ALGORITHM_FACTORY( QuarticPhotonReco )
 
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h
index 7f5d5fe2a010003617ada530115cbc3fbc371a40..a5a3b1e54d3de1285bc3c6fa9624e3dc8f4559f5 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h
@@ -36,11 +36,17 @@
 #include "RichDet/DeRichSphMirror.h"
 
 // Quartic Solver
-#include "RichRecUtils/QuarticSolver.h"
+#include "RichRecUtils/QuarticSolverNewton.h"
 
 // interfaces
 #include "RichInterfaces/IRichMirrorSegFinderLookUpTable.h"
 
+// Vector Class
+#include "VectorClass/instrset.h"
+
+// VDT
+#include "vdt/asin.h"
+
 namespace Rich
 {
   namespace Future
@@ -51,6 +57,9 @@ namespace Rich
       // Use the functional framework
       using namespace Gaudi::Functional;
 
+      // pull in methods from Rich::RayTracingUtils
+      using namespace Rich::RayTracingUtils;
+
       namespace
       {
         /// Shortcut to the output data type
@@ -98,8 +107,73 @@ namespace Rich
                             const Rich::PDPixelCluster::Vector& clusters,
                             const SpacePointVector& globalPoints,
                             const SpacePointVector& localPoints,
-                            const Relations::TrackToSegments::Vector& tkToSegRels ) const override;
+                            const Relations::TrackToSegments::Vector& tkToSegRels ) const override
+        {
+          return (this->*m_run)( segments, ckAngles, ckResolutions, trHitPntsLoc,
+                                 segPhotFlags, clusters, globalPoints, localPoints, tkToSegRels );
+        }
 
+      private: // dispatch methods
+
+        /** Pointer to dispatch function for the operator call
+         *  Default to lowest common denominator. 
+         *  Reset in initialize() according to runtime support */
+        OutData (QuarticPhotonReco::*m_run)
+          ( const LHCb::RichTrackSegment::Vector& segments,
+            const CherenkovAngles::Vector& ckAngles,
+            const CherenkovResolutions::Vector& ckResolutions,
+            const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+            const SegmentPhotonFlags::Vector& segPhotFlags,
+            const Rich::PDPixelCluster::Vector& clusters,
+            const SpacePointVector& globalPoints,
+            const SpacePointVector& localPoints,
+            const Relations::TrackToSegments::Vector& tkToSegRels ) const 
+          = &QuarticPhotonReco::run_generic;
+        
+        /// Generic method
+        OutData run_generic( const LHCb::RichTrackSegment::Vector& segments,
+                             const CherenkovAngles::Vector& ckAngles,
+                             const CherenkovResolutions::Vector& ckResolutions,
+                             const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                             const SegmentPhotonFlags::Vector& segPhotFlags,
+                             const Rich::PDPixelCluster::Vector& clusters,
+                             const SpacePointVector& globalPoints,
+                             const SpacePointVector& localPoints,
+                             const Relations::TrackToSegments::Vector& tkToSegRels ) const;
+
+        /// SSE4 method
+        OutData run_sse4( const LHCb::RichTrackSegment::Vector& segments,
+                          const CherenkovAngles::Vector& ckAngles,
+                          const CherenkovResolutions::Vector& ckResolutions,
+                          const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                          const SegmentPhotonFlags::Vector& segPhotFlags,
+                          const Rich::PDPixelCluster::Vector& clusters,
+                          const SpacePointVector& globalPoints,
+                          const SpacePointVector& localPoints,
+                          const Relations::TrackToSegments::Vector& tkToSegRels ) const;
+         
+        /// AVX method
+        OutData run_avx( const LHCb::RichTrackSegment::Vector& segments,
+                         const CherenkovAngles::Vector& ckAngles,
+                         const CherenkovResolutions::Vector& ckResolutions,
+                         const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                         const SegmentPhotonFlags::Vector& segPhotFlags,
+                         const Rich::PDPixelCluster::Vector& clusters,
+                         const SpacePointVector& globalPoints,
+                         const SpacePointVector& localPoints,
+                         const Relations::TrackToSegments::Vector& tkToSegRels ) const;
+
+        /// AVX2 method
+        OutData run_avx2( const LHCb::RichTrackSegment::Vector& segments,
+                          const CherenkovAngles::Vector& ckAngles,
+                          const CherenkovResolutions::Vector& ckResolutions,
+                          const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                          const SegmentPhotonFlags::Vector& segPhotFlags,
+                          const Rich::PDPixelCluster::Vector& clusters,
+                          const SpacePointVector& globalPoints,
+                          const SpacePointVector& localPoints,
+                          const Relations::TrackToSegments::Vector& tkToSegRels ) const;
+        
       private: // methods
 
         /// Access the DeRich beam pipe objects
@@ -115,15 +189,39 @@ namespace Rich
         inline float  myacos( const float  x ) const noexcept { return vdt::fast_acosf(x); }
 
         /// Find mirror segments and reflection points for given data
-        bool findMirrorData( const Rich::DetectorType rich,
-                             const Rich::Side side,
-                             const Gaudi::XYZPoint & virtDetPoint,
-                             const Gaudi::XYZPoint & emissionPoint,
-                             const DeRichSphMirror *& sphSegment,
-                             const DeRichSphMirror *& secSegment,
-                             Gaudi::XYZPoint & sphReflPoint,
-                             Gaudi::XYZPoint & secReflPoint ) const;
-
+        inline bool findMirrorData( const Rich::DetectorType rich,
+                                    const Rich::Side side,
+                                    const Gaudi::XYZPoint & virtDetPoint,
+                                    const Gaudi::XYZPoint & emissionPoint,
+                                    const DeRichSphMirror *& sphSegment,
+                                    const DeRichSphMirror *& secSegment,
+                                    Gaudi::XYZPoint & sphReflPoint,
+                                    Gaudi::XYZPoint & secReflPoint ) const
+        {
+          // solve quartic equation with nominal values and find spherical mirror reflection point
+          m_quarticSolver.solve<double,3,3>( emissionPoint,
+                                             m_rich[rich]->nominalCentreOfCurvature(side),
+                                             virtDetPoint,
+                                             m_rich[rich]->sphMirrorRadius(),
+                                             sphReflPoint );
+          // find the spherical mirror segment
+          sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
+          // Search for the secondary segment
+          // Direction vector from primary mirror point to virtual detection point
+          const auto dir ( virtDetPoint - sphReflPoint );
+          // find the sec mirror intersction point and secondary mirror segment
+          const bool OK = intersectPlane( sphReflPoint, dir,
+                                          m_rich[rich]->nominalPlane(side),
+                                          secReflPoint );
+          if ( OK )
+          {
+            // find the secondary mirror
+            secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint );
+          }
+          // return
+          return OK;
+        }
+        
         /// Compute the best emission point for the gas radiators using
         /// the given spherical mirror reflection points
         bool getBestGasEmissionPoint( const Rich::RadiatorType radiator,
@@ -136,8 +234,8 @@ namespace Rich
 
       private:
 
-        /// Quartic Solver
-        Rich::Rec::QuarticSolver m_quarticSolver;
+        /// Newton Quartic Solver 
+        Rich::Rec::QuarticSolverNewton m_quarticSolver;
 
         /// Rich1 and Rich2 detector elements
         DetectorArray<const DeRich *> m_rich = {{}};
@@ -224,8 +322,98 @@ namespace Rich
         ToolHandle<const IMirrorSegFinderLookUpTable> m_mirrorSegFinder
         { "Rich::Future::MirrorSegFinderLookUpTable/MirrorFinder:PUBLIC", this };
 
+        /// Flag to control if CPU detection should be performed or not
+        Gaudi::Property<bool> m_detectCPU { this, "DetectCPUCapabilites", true };
+
       };
 
+      //=========================================================================
+      // Compute the best emission point for the gas radiators using
+      // the given spherical mirror reflection points
+      //=========================================================================
+      inline bool
+      QuarticPhotonReco::
+      getBestGasEmissionPoint( const Rich::RadiatorType radiator,
+                               const Gaudi::XYZPoint& sphReflPoint1,
+                               const Gaudi::XYZPoint& sphReflPoint2,
+                               const Gaudi::XYZPoint& detectionPoint,
+                               const LHCb::RichTrackSegment & segment,
+                               Gaudi::XYZPoint & emissionPoint,
+                               float & fraction ) const
+      {
+        double alongTkFrac = 0.5;
+        
+        if ( radiator == Rich::Rich1Gas )
+        {
+          // First reflection and hit point on same y side ?
+          const bool sameSide1 = ( sphReflPoint1.y() * detectionPoint.y() > 0 );
+          const bool sameSide2 = ( sphReflPoint2.y() * detectionPoint.y() > 0 );
+          if ( sameSide1 && sameSide2 )
+          {
+            emissionPoint = segment.bestPoint();
+          }
+          else if ( sameSide1 )
+          {
+            fraction = (float)(std::fabs(sphReflPoint1.y()/(sphReflPoint1.y()-sphReflPoint2.y())));
+            alongTkFrac = fraction/2.0;
+            emissionPoint = segment.bestPoint(alongTkFrac);
+          }
+          else if ( sameSide2 )
+          {
+            fraction = (float)(std::fabs(sphReflPoint2.y()/(sphReflPoint1.y()-sphReflPoint2.y())));
+            alongTkFrac = 1.0-fraction/2.0;
+            emissionPoint = segment.bestPoint(alongTkFrac);
+          }
+          else
+          {
+            //Warning( "Rich1Gas : Both RICH spherical mirror hits opposite side to detection point" );
+            return false;
+          }
+          
+        }
+        else if ( radiator == Rich::Rich2Gas )
+        {
+          // First sphReflPoint and hit point on same x side ?
+          const bool sameSide1 = ( sphReflPoint1.x() * detectionPoint.x() > 0 );
+          const bool sameSide2 = ( sphReflPoint2.x() * detectionPoint.x() > 0 );
+          if ( sameSide1 && sameSide2 )
+          {
+            emissionPoint = segment.bestPoint();
+          }
+          else if ( sameSide1 )
+          {
+            fraction = (float)(std::fabs(sphReflPoint1.x()/(sphReflPoint1.x()-sphReflPoint2.x())));
+            alongTkFrac = fraction/2.0;
+            emissionPoint = segment.bestPoint(alongTkFrac);
+          }
+          else if ( sameSide2 )
+          {
+            fraction = (float)(std::fabs(sphReflPoint2.x()/(sphReflPoint1.x()-sphReflPoint2.x())));
+            alongTkFrac = 1.0-fraction/2.0;
+            emissionPoint = segment.bestPoint(alongTkFrac);
+          }
+          else
+          {
+            //Warning( "Rich2Gas : Both RICH spherical mirror hits opposite side to detection point" );
+            return false;
+          }
+        }
+        else { Error( "::getBestGasEmissionPoint() called for Aerogel segment !!" ); }
+        
+        // if ( msgLevel(MSG::VERBOSE) )
+        // {
+        //   verbose() << radiator << " best emission point correction :- " << endmsg
+        //             << " -> Photon detection point = " << detectionPoint << endmsg
+        //             << " -> Sph. Mirror ptns       = " << sphReflPoint1 << " " << sphReflPoint2 << endmsg
+        //             << " -> Segment entry/exit     = " << trSeg.entryPoint() << " " << trSeg.exitPoint() << endmsg
+        //             << " -> Segment fraction       = " << fraction << endmsg
+        //             << " -> Emm. Ptn. Along traj   = " << alongTkFrac << endmsg
+        //             << " -> Best Emission point    = " << emissionPoint << endmsg;
+        // }
+        
+        return true;
+      }
+
     }
   }
 }
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.icpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.icpp
new file mode 100644
index 0000000000000000000000000000000000000000..63eb84b8e528f4df75650fbd2bd017b21506af6e
--- /dev/null
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.icpp
@@ -0,0 +1,534 @@
+
+{
+  // use ranges v3 library
+  using namespace ranges::v3;
+
+  // make the output data
+  OutData outData;
+
+  // Shortcut to the photons and relations
+  auto & photons   = std::get<0>(outData);
+  auto & relations = std::get<1>(outData);
+  // guess at reserve size
+  const auto resSize = segments.size() * globalPoints.size() / 20;
+  photons.reserve( resSize );
+  relations.reserve( resSize );
+
+  // View range for pixel data
+  const auto allPixels = Ranges::ConstZip( clusters, globalPoints, localPoints );
+
+  // Make 'ranges' for RICH pixels
+  const auto RPixs = richRanges( allPixels );
+  const auto & RichRanges  = RPixs[0];   // Complete RICH ranges
+  const auto & Rich1Ranges = RPixs[1];   // ranges for RICH1 [top,bottom]
+  const auto & Rich2Ranges = RPixs[2];   // ranges for RICH2 [left,right]
+  
+  // local position corrector
+  // longer term need to remove this
+  const Rich::Rec::RadPositionCorrector corrector;
+
+  // global photon index
+  int photonIndex(-1);
+  
+  // Loop over the track->segment relations
+  for ( const auto & inTkRel : tkToSegRels )
+  {
+    // loop over segments for this track
+    for ( const auto & segIndex : inTkRel.segmentIndices )
+    {
+      // Get the data from the segment containers
+      const auto & segment     = segments[segIndex];
+      const auto & tkCkAngles  = ckAngles[segIndex];
+      const auto & tkCkRes     = ckResolutions[segIndex];
+      const auto & tkLocPtn    = trHitPntsLoc[segIndex];
+      const auto & segFlags    = segPhotFlags[segIndex];
+
+      // Is this segment above threshold. Should be implicit due to above segment loop
+      //if ( ! ( tkCkAngles[lightestActiveHypo()] > 0 ) ) continue;
+
+      // which RICH and radiator
+      const auto rich = segment.rich();
+      const auto rad  = segment.radiator();
+
+      // This implementiation does not support Aerogel
+      if ( UNLIKELY( Rich::Aerogel == rad ) ) { Exception("Aerogel not supported"); }
+
+      // get the best pixel range for this segment, based on where hits are expected
+      const auto& pixR = ( segFlags.inBothPanels() ? RichRanges[rich] :
+                           Rich::Rich1 == rich ? 
+                           ( segFlags.inPanel(Rich::top) ? 
+                             Rich1Ranges[Rich::top]  : Rich1Ranges[Rich::bottom] ) :
+                           ( segFlags.inPanel(Rich::left) ? 
+                             Rich2Ranges[Rich::left] : Rich2Ranges[Rich::right] ) );
+
+      // Loop over pixel information
+      int pixIndex = pixR.first - allPixels.begin() - 1; // (start index - 1) for this range
+      for ( auto pix_it = pixR.first; pix_it != pixR.second; ++pix_it )
+      {
+        // load the data from the tuple
+        const auto & cluster = std::get<0>(*pix_it);
+        const auto & gloPos  = std::get<1>(*pix_it);
+        const auto & locPos  = std::get<2>(*pix_it);
+        // increment the pixel index
+        ++pixIndex;
+
+        //_ri_verbo << "Trying Segment P=" << segment.bestMomentum() << " Pixel " 
+        //          << cluster << endmsg;
+
+        // Apply pre-selection
+      
+        // Pixel position, in local HPD coords corrected for average radiator distortion
+        // Might need to pre-compute and cache this as we repeat it here each segment...
+        const auto pixP = corrector.correct(locPos,rad);
+
+        //  Track local hit point on the same panel as the hit
+        /// @todo CRJ Need to figure out why I was not correcting this in the current stack...
+        ///           Probably correct but should check. 
+        //const auto segPanelPnt = corrector.correct(tkLocPtn.point(cluster.panel()),rad);
+        const auto& segPanelPnt = tkLocPtn.point(cluster.panel());
+
+        // compute the seperation squared
+        const auto sep2 = ( std::pow( (pixP.x()-segPanelPnt.x()), 2 ) +
+                            std::pow( (pixP.y()-segPanelPnt.y()), 2 ) );
+
+        // presel. default to rejected
+        bool SelOK = false;
+
+        // Check overall boundaries
+        //_ri_verbo << " -> sep2 = " << sep2 << " " << m_maxROI2PreSel[rad]
+        //          << " " << m_minROI2PreSel[rad] << endmsg;
+        if ( sep2 < m_maxROI2PreSel[rad] && sep2 > m_minROI2PreSel[rad] )
+        {
+        
+          // estimated CK theta
+          const auto ckThetaEsti = std::sqrt(sep2)*m_scalePreSel[rad];
+          //_ri_verbo << " -> CK theta Esti " << ckThetaEsti << endmsg;
+
+          // Is the hit close to any mass hypo in local coordinate space
+          SelOK = any_of( activeParticles(),
+                         [sigma=m_nSigmaPreSel[rad],&ckThetaEsti,&tkCkAngles,&tkCkRes]
+                         ( const auto hypo )
+                         { return fabs(tkCkAngles[hypo]-ckThetaEsti) < sigma*tkCkRes[hypo]; } );
+          
+        } // boundary check
+
+        // Did we pass the pre-sel
+        if ( !SelOK )
+        {
+          //_ri_verbo << "   -> FAILED pre-selection -> reject" << endmsg;
+          continue;
+        }
+
+        // Move on to the full Quartic reconstruction
+
+        //std::cout << "Starting photon reco" << std::endl;
+
+        // make a photon object to work on
+        photons.emplace_back();
+        auto & gPhoton = photons.back();
+        DataGuard<decltype(photons)> photGuard(photons);
+
+        // pixel side
+        const auto side = cluster.panel();
+
+        // Emission point to use for photon reconstruction
+        // operate directly on photon data member for efficiency
+        //auto & emissionPoint = gPhoton.emissionPoint();
+        //emissionPoint = segment.bestPoint();
+        auto emissionPoint = segment.bestPoint();
+
+        // Photon direction at emission point
+        // Again, operator directly on data member
+        //auto & photonDirection = gPhoton.emissionDir();
+        Gaudi::XYZVector photonDirection;
+
+        // Final reflection points on sec and spherical mirrors
+        // operate directly on photon data
+        //auto & sphReflPoint = gPhoton.sphMirReflectionPoint();
+        //auto & secReflPoint = gPhoton.flatMirReflectionPoint();
+        Gaudi::XYZPoint sphReflPoint, secReflPoint;
+
+        // fraction of segment path length accessible to the photon
+        float fraction(1);
+
+        // Pointers to best sec and spherical mirror segments
+        const DeRichSphMirror * sphSegment = nullptr;
+        const DeRichSphMirror * secSegment = nullptr;
+
+        // flag to say if this photon candidate is un-ambiguous - default to false
+        bool unambigPhoton = false;
+
+        // find the reflection of the detection point in the sec mirror
+        // (virtual detection point) starting with nominal values
+        // At this we are assuming a flat nominal mirror common to all segments
+        const auto nom_dist = m_rich[rich]->nominalPlane(side).Distance(gloPos);
+        auto virtDetPoint = gloPos - 2.0 * nom_dist * m_rich[rich]->nominalNormal(side);
+
+        // --------------------------------------------------------------------------------------
+        // For gas radiators, try start and end points to see if photon is unambiguous
+        // NOTE : For this we are using the virtual detection point determined using
+        // the noimnal flat secondary mirror plane. Now the secondary mirrors are actually
+        // spherical this may be introducing some additional uncertainties.
+        // --------------------------------------------------------------------------------------
+        if ( m_testForUnambigPhots[rad] )
+        {
+          if ( UNLIKELY( rad == Rich::Aerogel ) )
+          {
+            // use default emission point and assume unambiguous since path length is so short..
+            unambigPhoton = true;
+          }
+          else  // gas radiators
+          {
+          
+            // flag for beam pipe check. Default is OK
+            bool beamTestOK(true);
+          
+            // -------------------------------------------------------------------------------
+            // First emission point, at start of track segment
+            const auto emissionPoint1 = segment.bestPoint(0.01);
+            // Find mirror segments for this emission point
+            const DeRichSphMirror* sphSegment1 = nullptr;
+            const DeRichSphMirror* secSegment1 = nullptr;
+            Gaudi::XYZPoint sphReflPoint1, secReflPoint1;
+            if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint1,
+                                  sphSegment1, secSegment1, sphReflPoint1, secReflPoint1 ) )
+            {
+              //_ri_debug << rad << " : Failed to reconstruct photon for start of segment" << endmsg;
+              continue;
+            }
+            if ( m_checkBeamPipe[rad] )
+            {
+              beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint1, sphReflPoint1 );
+            }
+            // -------------------------------------------------------------------------------
+          
+            // -------------------------------------------------------------------------------
+            // now do it again for emission point #2, at end of segment
+            const auto emissionPoint2 = segment.bestPoint(0.99);
+            // Find mirror segments for this emission point
+            const DeRichSphMirror* sphSegment2 = nullptr;
+            const DeRichSphMirror* secSegment2 = nullptr;
+            Gaudi::XYZPoint sphReflPoint2, secReflPoint2;
+            if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint2,
+                                  sphSegment2, secSegment2, sphReflPoint2, secReflPoint2 ) )
+            {
+              //_ri_debug << rad << " : Failed to reconstruct photon for end of segment" << endmsg;
+              continue;
+            }
+            if ( !beamTestOK && m_checkBeamPipe[rad] )
+            {
+              beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint2, sphReflPoint2 );
+            }
+            // -------------------------------------------------------------------------------
+          
+            // -------------------------------------------------------------------------------
+            if ( !beamTestOK )
+            {
+              // both start and end points failed beam pipe test -> reject
+              //_ri_debug << rad << " : Failed ambiguous photon beampipe intersection checks" << endmsg;
+              continue;
+            }
+            // -------------------------------------------------------------------------------
+          
+            // -------------------------------------------------------------------------------
+            // Test to see if photon direction hits the real physical mirror segments
+            // -------------------------------------------------------------------------------
+            if ( UNLIKELY(m_checkPrimMirrSegs[rad]) )
+            {
+              // primary mirrors
+              const auto ok = ( sphSegment1->intersects( emissionPoint1,
+                                                         sphReflPoint1-emissionPoint1 ) ||
+                                sphSegment2->intersects( emissionPoint2,
+                                                         sphReflPoint2-emissionPoint2 ) );
+              if ( !ok )
+              {
+                //_ri_debug << rad << " : Failed mirror segment intersection checks" << endmsg;
+                continue;
+              }
+            }
+            // -------------------------------------------------------------------------------
+          
+            // -------------------------------------------------------------------------------
+            // Get the best gas emission point
+            if ( !getBestGasEmissionPoint( rad,
+                                           sphReflPoint1,
+                                           sphReflPoint2,
+                                           gloPos,
+                                           segment,
+                                           emissionPoint,
+                                           fraction ) )
+            {
+              //_ri_debug << rad << " : Failed to compute best gas emission point" << endmsg;
+              continue;
+            }
+            // -------------------------------------------------------------------------------
+          
+            // -------------------------------------------------------------------------------
+            // Is this an unambiguous photon - I.e. only one possible mirror combination
+            if ( ( sphSegment1 == sphSegment2 ) && ( secSegment1 == secSegment2 ) )
+            {
+              // Set pointers to the mirror detector objects
+              sphSegment = sphSegment1;
+              secSegment = secSegment1;
+              // rough guesses at reflection points (improved later on)
+              sphReflPoint = sphReflPoint1 + 0.5*(sphReflPoint2-sphReflPoint1);
+              secReflPoint = secReflPoint1 + 0.5*(secReflPoint2-secReflPoint1);
+              // photon is not unambiguous
+              unambigPhoton = true;
+            }
+            // -------------------------------------------------------------------------------
+          
+          } // end radiator type if
+        
+          // if configured to do so reject ambiguous photons
+          if ( UNLIKELY( m_rejectAmbigPhots[rad] && !unambigPhoton ) )
+          {
+            //_ri_debug << rad << " : Failed ambiguous photon test" << endmsg;
+            continue;
+          }
+        
+        } // end unambiguous photon check
+        // --------------------------------------------------------------------------------------
+
+        // --------------------------------------------------------------------------------------
+        // Active segment fraction cut
+        // --------------------------------------------------------------------------------------
+        if ( UNLIKELY( fraction < m_minActiveFrac[rad] ) )
+        {
+          //_ri_debug << rad << " : Failed active segment fraction cut" << endmsg;
+          continue;
+        }
+        // --------------------------------------------------------------------------------------
+
+        // --------------------------------------------------------------------------------------
+        // if ambiguous gas photon or if ambiguous photon check above has been skipped, try again
+        // using best emission point and nominal mirror geometries to get the spherical and sec
+        // mirrors. Also, force this reconstruction if the above unambiguous test was skipped
+        // --------------------------------------------------------------------------------------
+        if ( UNLIKELY( !m_testForUnambigPhots[rad] || !unambigPhoton ) )
+        {
+          if ( !findMirrorData(rich,side,
+                               virtDetPoint,emissionPoint,
+                               sphSegment,secSegment,sphReflPoint,secReflPoint ) )
+          {
+            //_ri_debug << rad << " : Failed backup photon reconstruction" << endmsg;
+            continue;
+          }
+        }
+        // --------------------------------------------------------------------------------------
+      
+        // --------------------------------------------------------------------------------------
+        // Finally reconstruct the photon using best emission point and the best mirror segments
+        // --------------------------------------------------------------------------------------
+        if ( m_useAlignedMirrSegs[rad] )
+        {
+        
+          // If iterations are disabled for this radiator, update the virtual detection point
+          // using the selected mirror segment
+          if ( UNLIKELY( 0 == m_nMaxQits[rad] ) )
+          {
+          
+            // Form the virtual detection point
+            const auto distance = secSegment->centreNormalPlane().Distance(gloPos);
+            virtDetPoint = gloPos - 2.0 * distance * secSegment->centreNormal();
+          
+            // solve the quartic using the new data
+            m_quarticSolver.solve<double,3,4>( emissionPoint,
+                                               sphSegment->centreOfCurvature(),
+                                               virtDetPoint,
+                                               sphSegment->radius(),
+                                               sphReflPoint );
+          
+            // (re)find the spherical mirror segment
+            // CRJ - Is this needed ?
+            sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
+          
+          }
+          else
+          {
+          
+            // Iterate to final solution, improving the secondary mirror info
+            int iIt(0);
+            Gaudi::XYZPoint last_mirror_point(0,0,0);
+            bool loopOK = true;
+            while ( iIt < m_nMaxQits[rad] )
+            {
+            
+              // Get secondary mirror reflection point,
+              // using the best actual secondary mirror segment at this point
+              const auto dir = virtDetPoint - sphReflPoint;
+              const auto sc = ( m_treatSecMirrsFlat[rich] ? 
+                                intersectPlane( sphReflPoint, dir,
+                                                secSegment->centreNormalPlane(),
+                                                secReflPoint ) :
+                                intersectSpherical( sphReflPoint, dir,
+                                                    secSegment->centreOfCurvature(),
+                                                    secSegment->radius(),
+                                                    secReflPoint ) );
+              if ( UNLIKELY(!sc) )
+              {
+                //_ri_debug << rad << " : Failed to intersect nominal secondary mirror plane" << endmsg;
+                loopOK = false;
+                break;
+              }
+
+              // (re)find the secondary mirror
+              secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint );
+            
+              // Compute the virtual reflection point
+ 
+              // Construct plane tangential to secondary mirror passing through reflection point
+              //const Gaudi::Plane3D plane( secSegment->centreOfCurvature()-secReflPoint, secReflPoint );
+              // re-find the reflection of the detection point in the sec mirror
+              // (virtual detection point) with this mirror plane
+              //const auto distance     = plane.Distance(gloPos);
+              //virtDetPoint = gloPos - 2.0 * distance * plane.Normal();
+
+              // Same as above, just written out by hand and simplified a bit
+              const auto normV = secSegment->centreOfCurvature() - secReflPoint;
+              const auto distance = ( ( normV.X() * ( gloPos.X() - secReflPoint.X() ) ) + 
+                                      ( normV.Y() * ( gloPos.Y() - secReflPoint.Y() ) ) + 
+                                      ( normV.Z() * ( gloPos.Z() - secReflPoint.Z() ) ) );
+              virtDetPoint = gloPos - ( normV * ( 2.0 * distance / normV.Mag2() ) );
+
+              // solve the quartic using the new data
+              m_quarticSolver.solve<double,3,4>( emissionPoint,
+                                                 sphSegment->centreOfCurvature(),
+                                                 virtDetPoint,
+                                                 sphSegment->radius(),
+                                                 sphReflPoint );
+
+              // (re)find the spherical mirror segment
+              sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
+
+              // for iterations after the first, see if we are still moving
+              // If not, abort iterations early
+              if ( iIt > m_nMinQits[rad] &&
+                   (last_mirror_point-secReflPoint).Mag2() < m_minSphMirrTolIt[rad] ) { break; }
+
+              // store last sec mirror point
+              last_mirror_point = secReflPoint;
+
+              // increment iteration counter
+              ++iIt;
+            }
+            // if above while loop failed, abort this photon
+            if ( UNLIKELY(!loopOK) ) { continue; }
+          
+          }
+        
+        }
+        // --------------------------------------------------------------------------------------
+      
+        // --------------------------------------------------------------------------------------
+        // check that spherical mirror reflection point is on the same side as detection point
+        // and (if configured to do so) photon does not cross between detector sides
+        // --------------------------------------------------------------------------------------
+        if ( UNLIKELY( !sameSide(rad,sphReflPoint,virtDetPoint) ) )
+        {
+          //_ri_debug << rad << " : Reflection point on wrong side" << endmsg;
+          continue;
+        }
+        if ( UNLIKELY( m_checkPhotCrossSides[rad] &&
+                       !sameSide(rad,sphReflPoint,emissionPoint) ) )
+        {
+          //_ri_debug << rad << " : Photon crosses between detector sides" << endmsg;
+          continue;
+        }
+        //else { continue; } // uncomment to select ONLY crossing photons
+        // --------------------------------------------------------------------------------------
+      
+        // --------------------------------------------------------------------------------------
+        // For as radiators if ambiguous photon checks are disabled (since this is
+        // already done for these photons during those checks), check if the photon would have
+        // intersected with the beampipe
+        // --------------------------------------------------------------------------------------
+        if ( UNLIKELY( m_checkBeamPipe[rad] && !m_testForUnambigPhots[rad] ) )
+        {
+          if ( UNLIKELY( deBeam(rich)->testForIntersection( emissionPoint, sphReflPoint ) ) )
+          {
+            //_ri_debug << rad << " : Failed final beampipe intersection checks" << endmsg;
+            continue;
+          }
+        }
+        // --------------------------------------------------------------------------------------
+
+        // --------------------------------------------------------------------------------------
+        // If using aligned mirror segments, get the final sec mirror reflection
+        // point using the best mirror segments available at this point
+        // For RICH2, use the spherical nature of the scondary mirrors
+        // For RICH1, where they are much flatter, assume complete flatness
+        // --------------------------------------------------------------------------------------
+        if ( m_useAlignedMirrSegs[rad] )
+        {
+          const auto dir = virtDetPoint - sphReflPoint;
+          const auto sc = ( m_treatSecMirrsFlat[rich] ? 
+                            intersectPlane( sphReflPoint, dir,
+                                            secSegment->centreNormalPlane(),
+                                            secReflPoint ) :
+                            intersectSpherical( sphReflPoint, dir,
+                                                secSegment->centreOfCurvature(),
+                                                secSegment->radius(),
+                                                secReflPoint ) );
+          if ( !sc )
+          {
+            //_ri_debug << rad << " : Failed final secondary mirror plane intersection" << endmsg;
+            continue;
+          }
+        }
+        // --------------------------------------------------------------------------------------
+
+        // --------------------------------------------------------------------------------------
+        // calculate the cherenkov angles using the photon and track vectors
+        // --------------------------------------------------------------------------------------
+        photonDirection = (sphReflPoint-emissionPoint).Unit();
+        float thetaCerenkov(0), phiCerenkov(0);
+        segment.angleToDirection( photonDirection, thetaCerenkov, phiCerenkov );
+        // --------------------------------------------------------------------------------------
+
+        //---------------------------------------------------------------------------------------
+        // Apply fudge factor correction for small biases in CK theta
+        //---------------------------------------------------------------------------------------
+        thetaCerenkov += ckThetaCorrection(rad);
+        //---------------------------------------------------------------------------------------
+
+        // --------------------------------------------------------------------------------------
+        // Final checks on the Cherenkov angles
+        // --------------------------------------------------------------------------------------
+        if ( UNLIKELY( !checkAngles( rad, tkCkAngles, tkCkRes,
+                                     thetaCerenkov, phiCerenkov ) ) )
+        {
+          //_ri_verbo << "    -> photon FAILED checkAngleInRange test" << endmsg;
+          continue;
+        }
+        // --------------------------------------------------------------------------------------
+
+        // --------------------------------------------------------------------------------------
+        // Set (remaining) photon parameters
+        // --------------------------------------------------------------------------------------
+        gPhoton.setCherenkovTheta         ( thetaCerenkov       );
+        gPhoton.setCherenkovPhi           ( phiCerenkov         );
+        gPhoton.setActiveSegmentFraction  ( fraction            );
+        //gPhoton.setDetectionPoint         ( gloPos              );
+        gPhoton.setSmartID                ( cluster.primaryID() );
+        gPhoton.setUnambiguousPhoton      ( unambigPhoton       );
+        //gPhoton.setPrimaryMirror          ( sphSegment          );
+        //gPhoton.setSecondaryMirror        ( secSegment          );
+        //_ri_verbo << "Created photon " << gPhoton << endmsg;
+        // --------------------------------------------------------------------------------------
+     
+        // If we get here, keep the just made photon
+        photGuard.setOK();
+        // Save relations
+        relations.emplace_back( ++photonIndex, pixIndex, segIndex, inTkRel.tkIndex );
+        //info() << relations.back() << endmsg;
+
+      }
+
+    }
+    
+  }
+  
+  //_ri_debug << "Created " << photons.size() << " photons" << endmsg;
+  return outData;
+}
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/avx/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/avx/RichQuarticPhotonReco.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1840557ea2ddd9b2e1403f05c293205a77c93823
--- /dev/null
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/avx/RichQuarticPhotonReco.cpp
@@ -0,0 +1,30 @@
+
+// local
+#include "../RichQuarticPhotonReco.h"
+
+// All code is in general Rich reconstruction namespace
+using namespace Rich::Future;
+using namespace Rich::Future::Rec;
+
+struct avx_guard {
+  // see Agner Fog's optimization guide, 12.1 about mixing AVX and non-AVX code,
+  // (http://www.agner.org/optimize/optimizing_cpp.pdf)
+  // and preserving the YMM register state.
+  // Invoking __mm256_zeroupper seems to reduce the overhead when switching.
+  ~avx_guard() { _mm256_zeroupper(); }
+};
+
+OutData 
+QuarticPhotonReco::run_avx( const LHCb::RichTrackSegment::Vector& segments,
+                            const CherenkovAngles::Vector& ckAngles,
+                            const CherenkovResolutions::Vector& ckResolutions,
+                            const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                            const SegmentPhotonFlags::Vector& segPhotFlags,
+                            const Rich::PDPixelCluster::Vector& clusters,
+                            const SpacePointVector& globalPoints,
+                            const SpacePointVector& localPoints,
+                            const Relations::TrackToSegments::Vector& tkToSegRels ) const
+{
+  avx_guard guard{};
+#include "../RichQuarticPhotonReco.icpp"
+}
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/avx2/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/avx2/RichQuarticPhotonReco.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..672aa639912aeaa8181868a4cf4af04f680a54d5
--- /dev/null
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/avx2/RichQuarticPhotonReco.cpp
@@ -0,0 +1,30 @@
+
+// local
+#include "../RichQuarticPhotonReco.h"
+
+// All code is in general Rich reconstruction namespace
+using namespace Rich::Future;
+using namespace Rich::Future::Rec;
+
+struct avx_guard {
+  // see Agner Fog's optimization guide, 12.1 about mixing AVX and non-AVX code,
+  // (http://www.agner.org/optimize/optimizing_cpp.pdf)
+  // and preserving the YMM register state.
+  // Invoking __mm256_zeroupper seems to reduce the overhead when switching.
+  ~avx_guard() { _mm256_zeroupper(); }
+};
+
+OutData 
+QuarticPhotonReco::run_avx2( const LHCb::RichTrackSegment::Vector& segments,
+                             const CherenkovAngles::Vector& ckAngles,
+                             const CherenkovResolutions::Vector& ckResolutions,
+                             const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                             const SegmentPhotonFlags::Vector& segPhotFlags,
+                             const Rich::PDPixelCluster::Vector& clusters,
+                             const SpacePointVector& globalPoints,
+                             const SpacePointVector& localPoints,
+                             const Relations::TrackToSegments::Vector& tkToSegRels ) const
+{
+  avx_guard guard{};
+#include "../RichQuarticPhotonReco.icpp"
+}
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/generic/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/generic/RichQuarticPhotonReco.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d5db58c9144febf65e7adeec47cd82df50ae8939
--- /dev/null
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/generic/RichQuarticPhotonReco.cpp
@@ -0,0 +1,21 @@
+
+// local
+#include "../RichQuarticPhotonReco.h"
+
+// All code is in general Rich reconstruction namespace
+using namespace Rich::Future;
+using namespace Rich::Future::Rec;
+
+OutData 
+QuarticPhotonReco::run_generic( const LHCb::RichTrackSegment::Vector& segments,
+                                const CherenkovAngles::Vector& ckAngles,
+                                const CherenkovResolutions::Vector& ckResolutions,
+                                const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                                const SegmentPhotonFlags::Vector& segPhotFlags,
+                                const Rich::PDPixelCluster::Vector& clusters,
+                                const SpacePointVector& globalPoints,
+                                const SpacePointVector& localPoints,
+                                const Relations::TrackToSegments::Vector& tkToSegRels ) const
+{
+#include "../RichQuarticPhotonReco.icpp"
+}
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/sse4/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/sse4/RichQuarticPhotonReco.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..82131fc1330fee441861b78481e6720765de93bf
--- /dev/null
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/sse4/RichQuarticPhotonReco.cpp
@@ -0,0 +1,21 @@
+
+// local
+#include "../RichQuarticPhotonReco.h"
+
+// All code is in general Rich reconstruction namespace
+using namespace Rich::Future;
+using namespace Rich::Future::Rec;
+
+OutData 
+QuarticPhotonReco::run_sse4( const LHCb::RichTrackSegment::Vector& segments,
+                             const CherenkovAngles::Vector& ckAngles,
+                             const CherenkovResolutions::Vector& ckResolutions,
+                             const SegmentPanelSpacePoints::Vector& trHitPntsLoc,
+                             const SegmentPhotonFlags::Vector& segPhotFlags,
+                             const Rich::PDPixelCluster::Vector& clusters,
+                             const SpacePointVector& globalPoints,
+                             const SpacePointVector& localPoints,
+                             const Relations::TrackToSegments::Vector& tkToSegRels ) const
+{
+#include "../RichQuarticPhotonReco.icpp"
+}
diff --git a/Rich/RichFutureRecSys/examples/Brunel.py b/Rich/RichFutureRecSys/examples/Brunel.py
index a522b3d5058f1198d5b1860c74c74c1417ff9ae3..c6b69a4c3240d416ad4b05067c5500275a2511cb 100644
--- a/Rich/RichFutureRecSys/examples/Brunel.py
+++ b/Rich/RichFutureRecSys/examples/Brunel.py
@@ -45,15 +45,15 @@ Brunel().MCCheckSequence = richs+["PROTO"]
 #############################################################################
 
 # Timestamps in messages
-LHCbApp().TimeStamp = True
+#LHCbApp().TimeStamp = True
 
 Brunel().OutputType = 'None'
 importOptions("$APPCONFIGOPTS/Persistency/Compression-ZLIB-1.py")
 
-Brunel().EvtMax     = 5000
+Brunel().EvtMax     = 10000
 Brunel().PrintFreq  = 500
 
-#Brunel().Histograms = "Expert"
+Brunel().Histograms = "Expert"
 
 #ApplicationMgr().ExtSvc += [ "AuditorSvc" ]
 
diff --git a/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py b/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py
index 93f5018ebebef6b0bf6c1a834b1a45a11fd0f93e..c56359d4d38a201539cafcc9bb41cc0392b84854 100644
--- a/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py
+++ b/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py
@@ -21,3 +21,6 @@ from Configurables import Brunel, LHCbApp
 Brunel().DataType  = "2016"
 Brunel().Simulation = False
 Brunel().WithMC     = False
+
+#Brunel().DDDBtag   = "dddb-20150724"
+#Brunel().CondDBtag = "cond-20161011"
diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py
index 0186612021d870948756bb5d8a4cb2263f1e10a5..43d45a745e14b777f2824863308e2766ab7ce2ae 100644
--- a/Rich/RichFutureRecSys/examples/RichFuture.py
+++ b/Rich/RichFutureRecSys/examples/RichFuture.py
@@ -94,7 +94,7 @@ RichCheck = RichRecCheckers( inputTrackLocations = tkLocs,
 # The final sequence to run
 all = GaudiSequencer( "All", Members = [ fetcher, odinDecode, richDecode, 
                                          tkUnpack, tkFilt, RichRec ] )
-all.Members += [ RichMoni ]
+#all.Members += [ RichMoni ]
 #all.Members += [ mcPUnpack, RichCheck ]
 all.MeasureTime = True
 
diff --git a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py
index d4124f78e8cec8fe1fc426ab08a4d493ee41944d..c21660c283ff5cfc0226b1423bff3ed8e7fcff21 100644
--- a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py
+++ b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py
@@ -137,7 +137,7 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group.
                       preloadGeometry = False,
 
                       # Number points for mass hypothesis ring ray tracing
-                      NRingPoints = ( 100, 100, 100 ),
+                      NRingPoints = ( 96, 96, 96 ),
 
                       # Track Extrapolator type
                       trackExtrapolator = "TrackRungeKuttaExtrapolator",
@@ -243,9 +243,6 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group.
     if algprops["Detectors"][0] : pixname += "RICH1"
     if algprops["Detectors"][1] : pixname += "RICH2"
 
-    # The complete sequence
-    all = GaudiSequencer( "RichRecoSeq" + GroupName, MeasureTime = MeasureTime )
-
     # ==================================================================
     # RICH Pixel Treatment. Common to all track types and instanciations
     # of the RICH sequences.
@@ -291,8 +288,9 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group.
     pixels = GaudiSequencer( "RichPixels"+GroupName,
                              MeasureTime = MeasureTime,
                              Members = [ richCluster, richGloPoints, richLocPoints ] )
-    # Add to final sequence
-    all.Members += [ pixels ]
+        # The complete sequence
+    all = GaudiSequencer( "RichRecoSeq" + GroupName, MeasureTime = MeasureTime )
+    all.Members = [pixels]
 
     # Loop over tracks
     for tktype,trackLocation in inputTrackLocations.iteritems() :
@@ -308,6 +306,7 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group.
 
         # Sequencer for this track type
         tkSeq = GaudiSequencer( "Rich"+name+"Reco", MeasureTime = MeasureTime )
+        tkSeq.Members = [ ] # make sure empty
         all.Members += [ tkSeq ]
         
         # ==================================================================
diff --git a/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt b/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt
index d0eb77b7e8e02d7f9c925bea2cd98d97a0dab547..c4c1f41f2891f34a711319d20fec457c0315a3c2 100644
--- a/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt
+++ b/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt
@@ -15,10 +15,9 @@ find_package(Boost)
 find_package(GSL)
 find_package(VDT)
 
-include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS})
+include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} ${Vc_INCLUDE_DIR})
 
 gaudi_add_module(RichFutureRecTrackAlgorithms
                  src/*.cpp
                  INCLUDE_DIRS Boost GSL VDT Rich/RichFutureRecBase Rich/RichUtils Rich/RichFutureUtils Tr/TrackInterfaces Tr/TrackKernel
-                 LINK_LIBRARIES Boost GSL VDT RichFutureRecBase RichUtils RichFutureUtils)
-
+                 LINK_LIBRARIES Boost GSL VDT RichFutureRecBase RichUtils RichFutureUtils )
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp
index 14b9006693a24c30e7e2661373284d7326338442..65b1025f2b73c4393c208e6ce77a3a2affa4ab10 100755
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp
@@ -40,7 +40,7 @@ StatusCode BaseTrSegMaker::initialize()
 
   if ( !usedRads(Rich::Aerogel) )
   {
-    debug() << "Track segments for Aerogel are disabled" << endmsg;
+    _ri_debug << "Track segments for Aerogel are disabled" << endmsg;
     //Warning("Track segments for Aerogel are disabled",StatusCode::SUCCESS).ignore();
     //made into debug as it is still too verbose as a warning
   }
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h
index 2aa35f224f33551ac6af02b9eb4c83a1f515cbc7..f4dc7eaec57324037845f00f1ce7938400eb58a6 100755
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h
@@ -49,7 +49,7 @@ namespace Rich
     {
 
       //---------------------------------------------------------------------------------
-      /** @class BaseTrSegMakerFromRecoTracks RichBaseTrSegMakerFromRecoTracks.h
+      /** @class BaseTrSegMaker RichBaseTrSegMaker.h
        *
        *  Base class to tools that create RichTrackSegments from Tracks.
        *
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h
index 66423a3cf4a201d53b5b3b235a472c7bb596070e..67e050fba025351d76cc7855a5483e4243793740 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h
@@ -42,7 +42,7 @@ namespace Rich
                                     SegmentPhotonFlags::Vector >;
       }
 
-      /** @class RichGeomEffCKMassRings RichRichGeomEffCKMassRings.h
+      /** @class GeomEffCKMassRings RichRichGeomEffCKMassRings.h
        *
        *  Computes the geometrical efficiencies for a set of track segments.
        *
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp
index b39a6afa73fec8fa9a164047f61a2b48c9500f78..64cc4a39e0ff7e59e01987f3cc502ec83dbcf6ca 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp
@@ -110,46 +110,106 @@ RayTraceCherenkovCones::operator()( const LHCb::RichTrackSegment::Vector& segmen
         // reserve size in the points container
         rings[id].reserve( m_nPoints[rad] );
 
-        // loop around the ring
-        unsigned int nOK(0), nPhot(0);
-        for ( auto iP = m_cosSinPhiV[rad].begin(); iP != m_cosSinPhiV[rad].end(); ++iP, ++nPhot )
+        // Container for photon directions
+        std::vector<Gaudi::XYZVector> vects;
+        vects.reserve( m_nPoints[rad] );
+
+        // loop around the ring to create the directions
+        for ( const auto& P : m_cosSinPhiV[rad] )
         {
           // Photon direction around loop
-          const auto photDir =
-            segment.vectorAtCosSinThetaPhi( cosTheta,     sinTheta,
-                                            (*iP).cosPhi, (*iP).sinPhi );
-
-          // do the ray tracing
-          // really need to speed this up by removing virtual function call and need for photon. 
-          const auto result =
-            m_rayTrace.get()->traceToDetector( rich, emissionPoint, photDir, photon, 
+          vects.emplace_back( segment.vectorAtCosSinThetaPhi( cosTheta, sinTheta,
+                                                              P.cosPhi, P.sinPhi ) );
+        }
+
+        // Count the number of good photons
+        unsigned int nOK(0);
+
+        // Which ray tracing to run
+        if ( m_useVectorised )
+        {
+
+          // The vectorised ray tracing
+          const auto results =
+            m_rayTrace.get()->traceToDetector( emissionPoint, vects, segment, m_traceModeRad[rad] );
+
+          // loop over the results and fill
+          unsigned int nPhot(0);
+          for ( const auto && data : Ranges::ConstZip(results,m_cosSinPhiV[rad]) )
+          {
+            const auto & res    = std::get<0>(data);
+            const auto & cosphi = std::get<1>(data);
+
+            // count photons
+            ++nPhot;
+
+            // Add a new point
+            const auto & gP = res.detectionPoint;
+            rings[id].emplace_back ( gP,
+                                     corrector.correct(m_idTool.get()->globalToPDPanel(gP),rad),
+                                     res.smartID,
+                                     (RayTracedCKRingPoint::Acceptance)(res.result),
+                                     res.primaryMirror,
+                                     res.secondaryMirror,
+                                     cosphi.phi );
+
+            // count raytraces that are in HPD panel
+            if ( res.result >= LHCb::RichTraceMode::InHPDPanel ) { ++nOK; }
+
+            // bailout check
+            if ( 0 == nOK && nPhot >= m_nBailout[rad] ) break; 
+            
+          }
+
+        }
+        else
+        {
+
+          // Loop over the directions and ray trace them using the scalar method
+          unsigned int nPhot(0);
+          for ( const auto && data : Ranges::ConstZip(vects,m_cosSinPhiV[rad]) )
+          {
+            const auto & photDir = std::get<0>(data);
+            const auto & cosphi  = std::get<1>(data);
+
+            // count photons
+            ++nPhot;
+
+            // do the ray tracing
+            // really need to speed this up by removing virtual function call and need for photon. 
+            const auto result =
+              m_rayTrace.get()->traceToDetector( rich, emissionPoint, photDir, photon, 
                                                segment, 
                                                m_traceModeRad[rad], Rich::top ); // Note forced side is not used..
 
-          // Add a new point
-          const auto & gP = photon.detectionPoint();
-          rings[id].emplace_back ( gP,
-                                   corrector.correct(m_idTool.get()->globalToPDPanel(gP),rad),
-                                   photon.smartID(),
-                                   (RayTracedCKRingPoint::Acceptance)(result),
-                                   photon.primaryMirror(),
-                                   photon.secondaryMirror(),
-                                   (*iP).phi );
-
-          // count raytraces that are in HPD panel
-          if ( result >= LHCb::RichTraceMode::InHPDPanel ) { ++nOK; }
-          
-          // bailout check
-          if ( 0 == nOK && nPhot >= m_nBailout[rad] ) break;
-          
+            // Add a new point
+            const auto & gP = photon.detectionPoint();
+            rings[id].emplace_back ( gP,
+                                     corrector.correct(m_idTool.get()->globalToPDPanel(gP),rad),
+                                     photon.smartID(),
+                                     (RayTracedCKRingPoint::Acceptance)(result),
+                                     photon.primaryMirror(),
+                                     photon.secondaryMirror(),
+                                     cosphi.phi );
+            
+            // count raytraces that are in HPD panel
+            if ( result >= LHCb::RichTraceMode::InHPDPanel ) { ++nOK; }
+            
+            // bailout check
+            if ( 0 == nOK && nPhot >= m_nBailout[rad] ) break; 
+          }
+
         }
 
+        // if no good hits empty the container
+        //if ( 0 == nOK ) { rings[id].clear(); } 
+        
       }
-
+      
     }
-
+    
   }
-
+  
   return ringsV;
 }
 
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h
index cb615571789f928efd44f83d4560fb66ece4c2ee..26dd924732da68f56166a598f8a1d9925ae5990f 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h
@@ -100,7 +100,7 @@ namespace Rich
 
         /// Number of points to ray trace on each ring, for each radiator
         Gaudi::Property< RadiatorArray<unsigned int> > m_nPoints
-        { this, "NRingPoints", { 100u, 100u, 100u } };
+        { this, "NRingPoints", { 96u, 96u, 96u } };
 
         /// Bailout number
         RadiatorArray<unsigned int> m_nBailout = {{}};
@@ -122,6 +122,10 @@ namespace Rich
         Gaudi::Property< RadiatorArray<float> > m_bailoutFrac
         { this, "BailoutFraction", { 0.75, 0.75, 0.75 } };
 
+        /// Use the scalar or vector ray tracing
+        Gaudi::Property<bool> m_useVectorised
+        { this, "UseVectorisedRayTracing", true };
+
         /// Ray tracing tool
         ToolHandle<const IRayTracing> m_rayTrace
         { "Rich::Future::RayTracing/RayTracing", this };
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h
index c8fc19ba0fa986d5687978f700ed08e362e521fc..3072c9e1147d043888cc3b189f2e28feb9f0282c 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h
@@ -53,7 +53,7 @@ namespace Rich
     namespace Rec
     {
 
-      /** @class RichTrackFunctionalCherenkovResolutions RichRichTrackFunctionalCherenkovResolutions.h
+      /** @class TrackFunctionalCherenkovResolutions RichTrackFunctionalCherenkovResolutions.h
        *
        *  Computes the expected Cherenkov resolutions for the given track segments
        *
diff --git a/Rich/RichRecPhotonTools/CMakeLists.txt b/Rich/RichRecPhotonTools/CMakeLists.txt
index c551830b1f52218da53fd8590664a146f8ecf67d..231f8c58a3e2acc3727de49df1b393c585490394 100644
--- a/Rich/RichRecPhotonTools/CMakeLists.txt
+++ b/Rich/RichRecPhotonTools/CMakeLists.txt
@@ -24,9 +24,3 @@ gaudi_add_module(RichRecPhotonTools
                  INCLUDE_DIRS Eigen Boost GSL VectorClass src
                  LINK_LIBRARIES Boost GSL RichDetLib GaudiKernel LHCbKernel LHCbMathLib RichKernelLib RichRecBase)
 
-if(GAUDI_BUILD_TESTS)
-  gaudi_add_executable(RichPhotonRecoTest
-                       src/application/PhotonReco/*.cpp
-                       INCLUDE_DIRS ROOT Eigen GSL VectorClass src
-                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel ROOT)
-endif()
diff --git a/Rich/RichRecPhotonTools/src/application/PhotonReco/main.cpp b/Rich/RichRecPhotonTools/src/application/PhotonReco/main.cpp
deleted file mode 100644
index 20a2eb747de69fa4c1f156fe1aecec8d63ada915..0000000000000000000000000000000000000000
--- a/Rich/RichRecPhotonTools/src/application/PhotonReco/main.cpp
+++ /dev/null
@@ -1,81 +0,0 @@
-
-// Quartic Solver
-#include "RichRecUtils/QuarticSolver.h"
-
-// STL
-#include <random>
-#include <vector>
-#include <iostream>
-#include <string>
-#include <typeinfo>
-
-// Make an instance of the quartic solver
-Rich::Rec::QuarticSolver qSolver;
-
-Gaudi::XYZPoint sphReflPoint;
-
-class Data
-{
-public:
-  typedef std::vector<Data> Vector;
-public:
-  Gaudi::XYZPoint emissPnt;
-  Gaudi::XYZPoint centOfCurv;
-  Gaudi::XYZPoint virtDetPoint;
-  double           radius;
-public:
-  Data() 
-  {
-    // randomn generator
-    static std::default_random_engine gen;
-    // Distributions for each member
-    static std::uniform_real_distribution<double> r_emiss_x(-800,800), r_emiss_y(-600,600), r_emiss_z(10000,10500);
-    static std::uniform_real_distribution<double> r_coc_x(-3000,3000), r_coc_y(-20,20),     r_coc_z(3300,3400);
-    static std::uniform_real_distribution<double> r_vdp_x(-3000,3000), r_vdp_y(-200,200),   r_vdp_z(8100,8200);
-    static std::uniform_real_distribution<float>  r_rad(8500,8600);
-    // setup data
-    emissPnt     = Gaudi::XYZPoint( r_emiss_x(gen), r_emiss_y(gen), r_emiss_z(gen) );
-    centOfCurv   = Gaudi::XYZPoint( r_coc_x(gen),   r_coc_y(gen),   r_coc_z(gen)   );
-    virtDetPoint = Gaudi::XYZPoint( r_vdp_x(gen),   r_vdp_y(gen),   r_vdp_z(gen)   );
-    radius       = r_rad(gen);
-  }
-};
-
-template< class TYPE >
-void solve( const Data& data )
-{
-  qSolver.solve<TYPE>( data.emissPnt, 
-                       data.centOfCurv, 
-                       data.virtDetPoint,
-                       data.radius, 
-                       sphReflPoint );
-}
-
-template< class TYPE >
-void solve( const Data::Vector & dataV )
-{
-  std::cout << "Solving Quartic Equation for " 
-            << typeid(TYPE).name() 
-            << " Photons ..." << std::endl;
-  // iterate over the data and solve it...
-  for ( const auto& data : dataV ) { solve<TYPE>(data); }
-}
-
-int main ( int /*argc*/, char** /*argv*/ )
-{
-  const unsigned int nPhotons = 1e6;
-  
-  Data::Vector dataV;
-  dataV.reserve( nPhotons );
-  // Construct the data to work on
-  std::cout << "Creating " << nPhotons << " random photons ..." << std::endl;
-  for ( unsigned int i = 0; i < nPhotons; ++i ) { dataV.push_back( Data() ); }
-
-  // run the solver for floats
-  solve<float>( dataV );
-
-  // run the solver for doubles
-  solve<double>( dataV );
-
-  return 0;
-}
diff --git a/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h b/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h
index b1618bab6d96e3435962dfab59578025bbc03127..6331185cce1bc015c04b0f4840ae1eb8d5a998b8 100755
--- a/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h
+++ b/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h
@@ -56,7 +56,7 @@
 #include "gsl/gsl_poly.h"
 
 // Quartic Solver
-#include "RichRecUtils/QuarticSolver.h"
+#include "RichRecUtils/QuarticSolverNewton.h"
 
 namespace Rich
 {
@@ -175,7 +175,7 @@ namespace Rich
       mutable const ISnellsLawRefraction * m_snellsLaw = nullptr;
 
       /// Quartic Solver
-      QuarticSolver m_quarticSolver;
+      QuarticSolverNewton m_quarticSolver;
 
       /** @brief Flag to indicate if the unambiguous photon test should be performed
        *  for each radiator
diff --git a/Rich/RichRecTests/CMakeLists.txt b/Rich/RichRecTests/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5a43b8ec44fc9b7ec830b5ab2b67558f9359979b
--- /dev/null
+++ b/Rich/RichRecTests/CMakeLists.txt
@@ -0,0 +1,109 @@
+################################################################################
+# Package: RichRecTests
+################################################################################
+gaudi_subdir(RichRecTests v1r0)
+
+gaudi_depends_on_subdirs(Det/RichDet
+                         GaudiKernel
+                         Kernel/LHCbKernel
+                         Kernel/LHCbMath
+                         Kernel/VectorClass
+                         Kernel/TemplatedGenVector
+                         Kernel/Vc
+                         Rich/RichRecUtils)
+
+#find_package(vectorclass)
+find_package(Boost)
+find_package(GSL)
+find_package(Eigen)
+find_package(Vc)
+find_package(ROOT COMPONENTS MathCore GenVector )
+
+message(STATUS "Using Vc ${Vc_INCLUDE_DIR} ${Vc_LIB_DIR}")
+#message(STATUS "Vc DEFINITIONS ${Vc_DEFINITIONS}")
+#message(STATUS "Vc COMPILE flags ${Vc_COMPILE_FLAGS}")
+#message(STATUS "Vc ARCHITECTURE flags ${Vc_ARCHITECTURE_FLAGS}")
+
+include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS}
+                           ${EIGEN_INCLUDE_DIRS} ${Vc_INCLUDE_DIR})
+
+if(GAUDI_BUILD_TESTS)
+
+  gaudi_add_executable(RichPhotonRecoTestSSE4
+                       src/PhotonReco/main_sse4.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichPhotonRecoTestSSE4 -lrt )
+  set_property(SOURCE src/PhotonReco/main_sse4.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -msse4.2 " )  
+
+  gaudi_add_executable(RichPhotonRecoTestAVX
+                       src/PhotonReco/main_avx.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichPhotonRecoTestAVX -lrt )
+  set_property(SOURCE src/PhotonReco/main_avx.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx " )  
+
+  gaudi_add_executable(RichPhotonRecoTestAVX2
+                       src/PhotonReco/main_avx2.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichPhotonRecoTestAVX2 -lrt )
+  # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here
+  #set_property(SOURCE src/PhotonReco/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx2 " ) 
+  set_property(SOURCE src/PhotonReco/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx " ) 
+
+  gaudi_add_executable(RichPhotonRecoTestAVX2FMA
+                       src/PhotonReco/main_avx2fma.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichPhotonRecoTestAVX2FMA -lrt )
+  # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here
+  #set_property(SOURCE src/PhotonReco/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mfma -mavx2 " )  
+  set_property(SOURCE src/PhotonReco/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx " )  
+
+
+
+  gaudi_add_executable(RichRayTracingTestSSE4
+                       src/RayTracing/main_sse4.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichRayTracingTestSSE4 "-lrt ${Vc_LIB_DIR}/libVc.a" )
+  set_property(SOURCE src/RayTracing/main_sse4.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -msse4.2 " )
+
+  gaudi_add_executable(RichRayTracingTestAVX
+                       src/RayTracing/main_avx.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichRayTracingTestAVX "-lrt ${Vc_LIB_DIR}/libVc.a" )
+  set_property(SOURCE src/RayTracing/main_avx.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx " )
+
+  gaudi_add_executable(RichRayTracingTestAVX2
+                       src/RayTracing/main_avx2.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichRayTracingTestAVX2 "-lrt ${Vc_LIB_DIR}/libVc.a" )
+  # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here
+  #set_property(SOURCE src/RayTracing/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx2 " )
+  set_property(SOURCE src/RayTracing/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx " )
+
+  gaudi_add_executable(RichRayTracingTestAVX2FMA
+                       src/RayTracing/main_avx2fma.cpp
+                       INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src
+                       LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib )
+  target_link_libraries( RichRayTracingTestAVX2FMA "-lrt ${Vc_LIB_DIR}/libVc.a" )
+  # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here
+  #set_property(SOURCE src/RayTracing/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mfma -mavx2 " )
+  set_property(SOURCE src/RayTracing/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx " )
+
+  if(LCG_COMP STREQUAL "gcc" OR
+      (BINARY_TAG_COMP_NAME STREQUAL "gcc" AND BINARY_TAG_COMP_VERSION VERSION_LESS "5"))
+    set_property(SOURCE src/PhotonReco/main_avx.cpp     APPEND_STRING PROPERTY COMPILE_FLAGS  " -fabi-version=0 " )
+    set_property(SOURCE src/PhotonReco/main_avx2.cpp    APPEND_STRING PROPERTY COMPILE_FLAGS  " -fabi-version=0 " )
+    set_property(SOURCE src/PhotonReco/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS  " -fabi-version=0 " )
+    set_property(SOURCE src/RayTracing/main_avx.cpp     APPEND_STRING PROPERTY COMPILE_FLAGS  " -fabi-version=0 " )
+    set_property(SOURCE src/RayTracing/main_avx2.cpp    APPEND_STRING PROPERTY COMPILE_FLAGS  " -fabi-version=0 " )
+    set_property(SOURCE src/RayTracing/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS  " -fabi-version=0 " )
+  endif()
+
+
+endif()
diff --git a/Rich/RichRecTests/src/PhotonReco/main.icpp b/Rich/RichRecTests/src/PhotonReco/main.icpp
new file mode 100644
index 0000000000000000000000000000000000000000..dbc074555b940527d629f5e3882824aa66bbfd69
--- /dev/null
+++ b/Rich/RichRecTests/src/PhotonReco/main.icpp
@@ -0,0 +1,124 @@
+
+// Gaudi
+//#include "GaudiKernel/Transform3DTypes.h"
+#include "GaudiKernel/Point3DTypes.h"
+#include "GaudiKernel/Vector3DTypes.h"
+
+// Quartic Solver
+#include "RichRecUtils/QuarticSolverNewton.h"
+
+// STL
+#include <random>
+#include <vector>
+#include <iostream>
+#include <string>
+#include <typeinfo>
+
+// Make an instance of the quartic solver
+Rich::Rec::QuarticSolverNewton qSolver;
+
+ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<float>  > sphReflPointF;
+ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<double> > sphReflPointD;
+
+unsigned long long int time_diff( struct timespec *start, 
+                                  struct timespec *stop ) 
+{
+  if ((stop->tv_nsec - start->tv_nsec) < 0) {
+    return 1000000000ULL * (stop->tv_sec - start->tv_sec - 1) +
+      stop->tv_nsec - start->tv_nsec + 1000000000ULL;
+  } else {
+    return 1000000000ULL * (stop->tv_sec - start->tv_sec) +
+      stop->tv_nsec - start->tv_nsec;
+  }
+}
+
+template< typename TYPE >
+class Data
+{
+public:
+  typedef std::vector<Data> Vector;
+public:
+  ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<TYPE> > emissPnt;
+  ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<TYPE> > centOfCurv;
+  ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<TYPE> > virtDetPoint;
+  TYPE radius;
+public:
+  Data() 
+  {
+    // randomn generator
+    static std::default_random_engine gen;
+    // Distributions for each member
+    static std::uniform_real_distribution<TYPE> r_emiss_x(-800,800), r_emiss_y(-600,600), r_emiss_z(10000,10500);
+    static std::uniform_real_distribution<TYPE> r_coc_x(-3000,3000), r_coc_y(-20,20),     r_coc_z(3300,3400);
+    static std::uniform_real_distribution<TYPE> r_vdp_x(-3000,3000), r_vdp_y(-200,200),   r_vdp_z(8100,8200);
+    static std::uniform_real_distribution<TYPE> r_rad(8500,8600);
+    // setup data
+    emissPnt     = { r_emiss_x(gen), r_emiss_y(gen), r_emiss_z(gen) };
+    centOfCurv   = { r_coc_x(gen),   r_coc_y(gen),   r_coc_z(gen)   };
+    virtDetPoint = { r_vdp_x(gen),   r_vdp_y(gen),   r_vdp_z(gen)   };
+    radius       = r_rad(gen);
+  }
+};
+
+template< class TYPE, class POINT >
+inline void solve( const Data<TYPE>& data, POINT & sphReflPoint )
+{
+  qSolver.solve<TYPE,3,2>( data.emissPnt, 
+                           data.centOfCurv, 
+                           data.virtDetPoint,
+                           data.radius, 
+                           sphReflPoint );
+}
+
+template< class TYPE, class POINT  >
+unsigned long long int __attribute__((noinline)) 
+solve( const typename Data<TYPE>::Vector & dataV, POINT & sphReflPoint )
+{
+  unsigned long long int best_dur{ 99999999999999999 };
+
+  timespec start, end;
+
+  const unsigned int nTests = 1000;
+  for ( unsigned int i = 0; i < nTests; ++i )
+  {
+    clock_gettime(CLOCK_MONOTONIC_RAW,&start);
+    // iterate over the data and solve it...
+    for ( const auto& data : dataV ) { solve<TYPE>(data,sphReflPoint); }
+    clock_gettime(CLOCK_MONOTONIC_RAW,&end);
+    const auto duration = time_diff(&start,&end);
+    if ( duration < best_dur ) { best_dur = duration; }
+  }
+
+  return best_dur ;
+}
+
+int main ( int /*argc*/, char** /*argv*/ )
+{
+  const unsigned int nPhotons = 1e5;
+  
+  Data<double>::Vector dataVD;
+  Data<float>::Vector dataVF;
+  dataVD.reserve( nPhotons );
+  dataVF.reserve( nPhotons );
+  // Construct the data to work on
+  std::cout << "Creating " << nPhotons << " random photons ..." << std::endl;
+  for ( unsigned int i = 0; i < nPhotons; ++i ) 
+  { 
+    dataVD.emplace_back( );
+    dataVF.emplace_back( ); 
+  }
+
+  // run the solver for doubles
+  auto dTime = solve<double>( dataVD, sphReflPointD );
+  std::cout << "Double " << dTime << std::endl;
+
+  // run the solver for floats
+  auto fTime = solve<float>( dataVF, sphReflPointF );
+  std::cout << "Float  " << fTime << std::endl;
+
+  // make sure we are not optimized away
+  asm volatile ("" : "+x"(fTime));
+  asm volatile ("" : "+x"(dTime));
+
+  return 0;
+}
diff --git a/Rich/RichRecTests/src/PhotonReco/main_avx.cpp b/Rich/RichRecTests/src/PhotonReco/main_avx.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/PhotonReco/main_avx.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecTests/src/PhotonReco/main_avx2.cpp b/Rich/RichRecTests/src/PhotonReco/main_avx2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/PhotonReco/main_avx2.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecTests/src/PhotonReco/main_avx2fma.cpp b/Rich/RichRecTests/src/PhotonReco/main_avx2fma.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/PhotonReco/main_avx2fma.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecTests/src/PhotonReco/main_sse4.cpp b/Rich/RichRecTests/src/PhotonReco/main_sse4.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/PhotonReco/main_sse4.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecTests/src/RayTracing/main.icpp b/Rich/RichRecTests/src/RayTracing/main.icpp
new file mode 100644
index 0000000000000000000000000000000000000000..8fb4bbf6928ceeb611dc2373fc4bdb2a26c3195e
--- /dev/null
+++ b/Rich/RichRecTests/src/RayTracing/main.icpp
@@ -0,0 +1,270 @@
+
+// STL
+#include <random>
+#include <vector>
+#include <iostream>
+#include <string>
+#include <typeinfo>
+
+// LHCbMath
+#include "LHCbMath/EigenTypes.h"
+//#include "LHCbMath/VectorClassTypes.h"
+
+// LHCb TemplatedGenVector
+#include "TemplatedGenVector/DisplacementVector3D.h"
+#include "TemplatedGenVector/PositionVector3D.h"
+#include "TemplatedGenVector/Plane3D.h"
+
+// Rich Utils
+#include "RichUtils/RichRayTracingUtils.h"
+#include "RichUtils/ZipRange.h"
+
+#include <Vc/Vc>
+
+unsigned long long int time_diff( struct timespec *start, 
+                                  struct timespec *stop ) 
+{
+  if ((stop->tv_nsec - start->tv_nsec) < 0) {
+    return 1000000000ULL * (stop->tv_sec - start->tv_sec - 1) +
+      stop->tv_nsec - start->tv_nsec + 1000000000ULL;
+  } else {
+    return 1000000000ULL * (stop->tv_sec - start->tv_sec) +
+      stop->tv_nsec - start->tv_nsec;
+  }
+}
+
+// randomn generator
+static std::default_random_engine gen;
+// Distributions for each member
+static std::uniform_real_distribution<double> p_x(-800,800),  p_y(-600,600), p_z(10000,10500);
+static std::uniform_real_distribution<double> d_x(-0.2,0.2),  d_y(-0.1,0.1), d_z(0.95,0.99);
+static std::uniform_real_distribution<double> c_x(3100,3200), c_y(10,15),    c_z(3200,3300);
+static std::uniform_real_distribution<double> r_rad(8500,8600);
+static std::uniform_real_distribution<double> p0(-0.002,0.002), p1(-0.2,0.2), p2(0.97,0.99), p3(-1300,1300);
+
+template < typename POINT, typename VECTOR, typename PLANE, typename FTYPE >
+class Data
+{
+public:
+  typedef std::vector< Data, Vc::Allocator<Data> > Vector;
+public:
+  POINT position;
+  VECTOR direction;
+  POINT CoC;
+  PLANE plane;
+  FTYPE radius{0};
+public:
+  template< typename INDATA >
+  Data( const INDATA & ind )
+    : position(ind.position),
+      direction(ind.direction),
+      CoC(ind.CoC),
+      plane(ind.plane),
+      radius(ind.radius) { }
+  Data( )
+    : position  ( p_x(gen), p_y(gen), p_z(gen) ),
+      direction ( d_x(gen), d_y(gen), d_z(gen) ),
+      CoC       ( c_x(gen), c_y(gen), c_z(gen) ),
+      plane     ( p0(gen),  p1(gen),  p2(gen), p3(gen) ), 
+      radius    ( r_rad(gen) ) { }
+};
+
+template < typename DATA >
+inline void trace( DATA & dataV )
+{
+  using namespace Rich::RayTracingUtils;
+  for ( auto & data : dataV )
+  {
+    reflectSpherical( data.position, data.direction, data.CoC, data.radius );
+    reflectPlane    ( data.position, data.direction, data.plane );
+    //std::cout << data.position << std::endl;
+  }
+}
+
+
+template < typename DATA >
+unsigned long long int __attribute__((noinline)) rayTrace( const DATA & in_dataV )
+{
+  timespec start, end;
+
+  unsigned long long int best_dur{ 99999999999999999 };
+
+  const unsigned int nTests = 10000;
+  //std::cout << "Running " << nTests << " tests on " << in_dataV.size() << " data objects" << std::endl;
+  DATA dataV; // working data copy;
+  for ( unsigned int i = 0; i < nTests; ++i )
+  {
+    dataV = in_dataV; // copy data, as it will be modified
+    clock_gettime(CLOCK_MONOTONIC_RAW,&start);
+    trace( dataV );
+    clock_gettime(CLOCK_MONOTONIC_RAW,&end);
+    const auto duration = time_diff(&start,&end);
+    if ( duration < best_dur ) { best_dur = duration; }
+  }
+
+  return best_dur ;
+}
+
+int main ( int /*argc*/, char** /*argv*/ )
+{
+
+  const unsigned int nPhotons = 96 * 1e2;
+  //const unsigned int nPhotons = 8;
+  std::cout << "Creating " << nPhotons << " random photons ..." << std::endl;
+
+  // Templated GenVector
+  using LHCbXYZPointD  = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<double>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using LHCbXYZVectorD = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<double>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using LHCbXYZPointF  = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<float>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using LHCbXYZVectorF = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<float>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using LHCbPlane3DD   = LHCbROOT::Math::Impl::Plane3D<double>;
+  using LHCbPlane3DF   = LHCbROOT::Math::Impl::Plane3D<float>;
+
+  // Vc
+  using VcXYZPointD  = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<Vc::double_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using VcXYZVectorD = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<Vc::double_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using VcXYZPointF  = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<Vc::float_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using VcXYZVectorF = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<Vc::float_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>;
+  using VcPlane3DD   = LHCbROOT::Math::Impl::Plane3D<Vc::double_v>;
+  using VcPlane3DF   = LHCbROOT::Math::Impl::Plane3D<Vc::float_v>;
+
+  //if ( false )
+  {
+
+    // Data in 'ROOT' SMatrix double format
+    Data<LHCbXYZPointD,LHCbXYZVectorD,LHCbPlane3DD,LHCbXYZPointD::Scalar>::Vector root_dataV( nPhotons );
+    auto resROOT = rayTrace( root_dataV );
+    std::cout << "ROOT double        " << resROOT << std::endl;
+
+    // Same data in Eigen double format
+    Data<LHCb::Math::Eigen::XYZPoint<double>,LHCb::Math::Eigen::XYZVector<double>,LHCb::Math::Eigen::Plane3D<double>,double>::Vector eigen_dataV( nPhotons );
+    auto resEIGEN = rayTrace( eigen_dataV );
+    std::cout << "Eigen double       " << resEIGEN << " speedup " << (float)resROOT/(float)resEIGEN << std::endl;
+
+    // Same data in VectorClass double format
+    //Data<LHCb::Math::VectorClass::XYZPoint<double>,LHCb::Math::VectorClass::XYZVector<double>,double>::Vector vclass_dataV( nPhotons );
+    //auto resVClass = rayTrace( vclass_dataV );
+    //std::cout << "VectorClass double " << resVClass << " speedup " << (float)resROOT/(float)resVClass << std::endl;
+
+    // Vc float
+    Data<VcXYZPointD,VcXYZVectorD,VcPlane3DD,Vc::double_v>::Vector VC_dataV( nPhotons / Vc::double_v::Size );
+    auto resVC = rayTrace( VC_dataV );
+    std::cout << "Vc double          " << resVC << " speedup " << (float)resROOT/(float)resVC << std::endl;
+
+    // make sure we are not optimized away
+    asm volatile ("" : "+x"(resROOT));
+    asm volatile ("" : "+x"(resEIGEN));
+    //asm volatile ("" : "+x"(resVClass));
+    asm volatile ("" : "+x"(resVC));
+
+  }
+
+  //if ( false )
+  {
+
+   // Data in 'ROOT' SMatrix float format
+    Data<LHCbXYZPointF,LHCbXYZVectorF,LHCbPlane3DF,LHCbXYZPointF::Scalar>::Vector root_dataV( nPhotons );
+    auto resROOT = rayTrace( root_dataV );
+    std::cout << "ROOT float         " << resROOT << std::endl;
+
+    // Same data in Eigen float format
+    Data<LHCb::Math::Eigen::XYZPoint<float>,LHCb::Math::Eigen::XYZVector<float>,LHCb::Math::Eigen::Plane3D<float>,float>::Vector eigen_dataV( nPhotons );
+    auto resEIGEN = rayTrace( eigen_dataV );
+    std::cout << "Eigen float        " << resEIGEN << " speedup " << (float)resROOT/(float)resEIGEN << std::endl;
+
+    // Same data in VectorClass float format
+    //Data<LHCb::Math::VectorClass::XYZPoint<float>,LHCb::Math::VectorClass::XYZVector<float>,float>::Vector vclass_dataV( nPhotons );
+    //auto resVClass = rayTrace( vclass_dataV );
+    //std::cout << "VectorClass float  " << resVClass << " speedup " << (float)resROOT/(float)resVClass << std::endl;
+
+    // Vc float
+    Data<VcXYZPointF,VcXYZVectorF,VcPlane3DF,Vc::float_v>::Vector VC_dataV( nPhotons / Vc::float_v::Size );
+    auto resVC = rayTrace( VC_dataV );
+    std::cout << "Vc float           " << resVC << " speedup " << (float)resROOT/(float)resVC << std::endl;
+
+    // make sure we are not optimized away
+    asm volatile ("" : "+x"(resROOT));
+    asm volatile ("" : "+x"(resEIGEN));
+    //asm volatile ("" : "+x"(resVClass));
+    asm volatile ("" : "+x"(resVC));
+
+  }
+
+  // Basic checks
+  if ( false )
+  {
+
+    std::vector< Vc::float_v > a = { 1, 2, 3, 4, 5, 6, 7, 8 };
+    std::vector< Vc::float_v > b = { 1, 2, 3, 4, 5, 6, 7, 8 };
+    
+    std::cout << "Size " << a.size() << " " << Vc::float_v::Size << std::endl;
+    
+    for ( unsigned int i = 0; i < (a.size()/Vc::float_v::Size); ++i )
+    {
+      a[i] += b[i];
+      std::cout << "Filling " << i << " " << a[i] << " " << b[i] << std::endl;
+    }
+    
+    for ( unsigned int i = 0; i < a.size(); ++i )
+    {
+      std::cout << "Reading " << i << " " << a[i] << " " << b[i] << std::endl;
+    }
+
+  }
+
+
+  // VC example
+  if ( false )
+  {
+    
+    using Vc::float_v;
+    //! [includes]
+    
+    //! [memory allocation]
+
+    // allocate memory for our initial x and y coordinates. Note that you can also put it into a
+    // normal float C-array but that you then must ensure alignment to Vc::VectorAlignment!
+    Vc::Memory<float_v, 100> x_mem;
+    Vc::Memory<float_v, 100> y_mem;
+    Vc::Memory<float_v, 100> r_mem;
+    Vc::Memory<float_v, 100> phi_mem;
+    //! [memory allocation]
+
+    //! [random init]
+    // fill the memory with values from -1.f to 1.f
+    std::cout << "vectorsCount = " << x_mem.vectorsCount() << std::endl;
+    for (size_t i = 0; i < x_mem.vectorsCount(); ++i)
+    {
+      x_mem.vector(i) = float_v::Random() * 2.f - 1.f;
+      y_mem.vector(i) = float_v::Random() * 2.f - 1.f;
+    }
+    //! [random init]
+
+    //! [conversion]
+    // calculate the polar coordinates for all coordinates and overwrite the euclidian coordinates
+    // with the result
+    for (size_t i = 0; i < x_mem.vectorsCount(); ++i) 
+    {
+      const float_v x = x_mem.vector(i);
+      const float_v y = y_mem.vector(i);
+      
+      r_mem.vector(i) = Vc::sqrt(x * x + y * y);
+      float_v phi = Vc::atan2(y, x) * 57.295780181884765625f; // 180/pi
+      phi(phi < 0.f) += 360.f;
+      phi_mem.vector(i) = phi;
+    }
+    //! [conversion]
+    
+    //! [output]
+    // print the results
+    std::cout << "entriesCount = " << x_mem.entriesCount() << std::endl;
+    for (size_t i = 0; i < x_mem.entriesCount(); ++i) {
+      std::cout << std::setw(3) << i << ": ";
+      std::cout << std::setw(10) << x_mem[i] << ", " << std::setw(10) << y_mem[i] << " -> ";
+      std::cout << std::setw(10) << r_mem[i] << ", " << std::setw(10) << phi_mem[i] << '\n';
+    }
+    
+  }
+
+  return 0;
+}
diff --git a/Rich/RichRecTests/src/RayTracing/main_avx.cpp b/Rich/RichRecTests/src/RayTracing/main_avx.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/RayTracing/main_avx.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecTests/src/RayTracing/main_avx2.cpp b/Rich/RichRecTests/src/RayTracing/main_avx2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/RayTracing/main_avx2.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecTests/src/RayTracing/main_avx2fma.cpp b/Rich/RichRecTests/src/RayTracing/main_avx2fma.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/RayTracing/main_avx2fma.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecTests/src/RayTracing/main_sse4.cpp b/Rich/RichRecTests/src/RayTracing/main_sse4.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e
--- /dev/null
+++ b/Rich/RichRecTests/src/RayTracing/main_sse4.cpp
@@ -0,0 +1,2 @@
+
+#include "main.icpp"
diff --git a/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h b/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h
index 4fbabe905e080b2efc52b75a53859659f7167d8f..067801e2169127a263cd49a9eed6e7d105459503 100755
--- a/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h
+++ b/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h
@@ -7,8 +7,7 @@
  */
 //----------------------------------------------------------------------
 
-#ifndef RICHRECPHOTONTOOLS_QuarticSolver_H
-#define RICHRECPHOTONTOOLS_QuarticSolver_H 1
+#pragma once
 
 // Gaudi
 #include "GaudiKernel/Kernel.h"
@@ -52,100 +51,6 @@ namespace Rich
     class QuarticSolver final
     {
 
-    public:
-
-      // Use eigen types
-      typedef LHCb::Math::Eigen::XYZPoint  Point;   ///< Point type
-      typedef LHCb::Math::Eigen::XYZVector Vector;  ///< vector type
-
-    public:
-
-      /** Solves the characteristic quartic equation for the RICH optical system.
-       *
-       *  See note LHCB/98-040 RICH section 3 for more details
-       *
-       *  @param emissionPoint Assumed photon emission point on track
-       *  @param CoC           Spherical mirror centre of curvature
-       *  @param virtDetPoint  Virtual detection point
-       *  @param radius        Spherical mirror radius of curvature
-       *  @param sphReflPoint  The reconstructed reflection pont on the spherical mirror
-       *
-       *  @return boolean indicating status of the quartic solution
-       *  @retval true  Calculation was successful. sphReflPoint is valid.
-       *  @retval false Calculation failed. sphReflPoint is not valid.
-       */
-      template< class TYPE >
-      inline void solve( const Gaudi::XYZPoint& emissionPoint,
-                         const Gaudi::XYZPoint& CoC,
-                         const Gaudi::XYZPoint& virtDetPoint,
-                         const TYPE radius,
-                         Gaudi::XYZPoint& sphReflPoint ) const
-      {
-
-        // typedefs vectorised types
-        typedef Eigen::Matrix< TYPE , 3 , 1 > Eigen3Vector;
-        using Vec4x = 
-          typename std::conditional<std::is_same<TYPE,float>::value,Vec4f,Vec4d>::type;
-
-        // vector from mirror centre of curvature to assumed emission point
-        const Vector evec( emissionPoint - CoC );
-        const TYPE e2 = evec.dot(evec);
-
-        // vector from mirror centre of curvature to virtual detection point
-        const Vector dvec( virtDetPoint - CoC );
-        const TYPE d2 = dvec.dot(dvec);
-
-        // various quantities needed to create quartic equation
-        // see LHCB/98-040 section 3, equation 3
-        const auto ed2 = e2 * d2;
-        const TYPE cosgamma2 = ( ed2 > 0 ? std::pow(evec.dot(dvec),2)/ed2 : 1.0 );
-        // vectorise 4 square roots into 1
-        const auto tmp_sqrt = sqrt( Vec4x( e2, d2,
-                                           cosgamma2 < 1.0 ? 1.0-cosgamma2 : 0.0,
-                                           cosgamma2 ) );
-        const auto e         = tmp_sqrt[0];
-        const auto d         = tmp_sqrt[1];
-        const auto singamma  = tmp_sqrt[2];
-        const auto cosgamma  = tmp_sqrt[3];
-
-        const auto dx        = d * cosgamma;
-        const auto dy        = d * singamma;
-        const auto r2        = radius * radius;
-        const auto dy2       = dy * dy;
-        const auto edx       = e + dx;
-
-        // Fill array for quartic equation
-        const auto a0      =     4.0 * ed2;
-        const auto inv_a0  =   ( a0 > 0 ? 1.0 / a0 : std::numeric_limits<TYPE>::max() );
-        const auto dyrad2  =     2.0 * dy * radius;
-        const auto a1      = - ( 2.0 * dyrad2 * e2 ) * inv_a0;
-        const auto a2      =   ( (dy2 * r2) + ( edx * edx * r2 ) - a0 ) * inv_a0;
-        const auto a3      =   ( dyrad2 * e * (e-dx) ) * inv_a0;
-        const auto a4      =   ( ( e2 - r2 ) * dy2 ) * inv_a0;
-
-        // use simplified RICH version of quartic solver
-        const auto sinbeta = solve_quartic_RICH<TYPE>( a1, // a
-                                                       a2, // b
-                                                       a3, // c
-                                                       a4  // d
-                                                       );
-
-        // (normalised) normal vector to reflection plane
-        auto n = evec.cross3(dvec);
-        n /= std::sqrt( n.dot(n) );
-
-        // construct rotation transformation
-        // Set vector magnitude to radius
-        // rotate vector and update reflection point
-        const Eigen::AngleAxis<TYPE> angleaxis( vdt::fast_asinf(sinbeta),
-                                                Eigen3Vector(n[0],n[1],n[2]) );
-        sphReflPoint = ( CoC + 
-                         Gaudi::XYZVector( angleaxis *
-                                           Eigen3Vector(evec[0],evec[1],evec[2]) *
-                                           ( radius / e ) ) );
-
-      }
-
     private:
 
       /// The cube root implementaton to use
@@ -169,10 +74,10 @@ namespace Rich
        */
       //----------------------------------------------------------------------
       template < class TYPE >
-      inline TYPE solve_quartic_RICH ( const TYPE& a,
-                                       const TYPE& b,
-                                       const TYPE& c,
-                                       const TYPE& d ) const 
+      inline TYPE solve_quartic_RICH ( const TYPE a,
+                                       const TYPE b,
+                                       const TYPE c,
+                                       const TYPE d ) const 
       {
 
         constexpr const TYPE r4 = 1.0 / 4.0;
@@ -237,9 +142,97 @@ namespace Rich
                  res );
       }
 
+    public:
+
+      /** Solves the characteristic quartic equation for the RICH optical system.
+       *
+       *  See note LHCB/98-040 RICH section 3 for more details
+       *
+       *  @param emissionPoint Assumed photon emission point on track
+       *  @param CoC           Spherical mirror centre of curvature
+       *  @param virtDetPoint  Virtual detection point
+       *  @param radius        Spherical mirror radius of curvature
+       *  @param sphReflPoint  The reconstructed reflection pont on the spherical mirror
+       *
+       *  @return boolean indicating status of the quartic solution
+       *  @retval true  Calculation was successful. sphReflPoint is valid.
+       *  @retval false Calculation failed. sphReflPoint is not valid.
+       */
+      template< class TYPE >
+      inline void solve( const Gaudi::XYZPoint& emissionPoint,
+                         const Gaudi::XYZPoint& CoC,
+                         const Gaudi::XYZPoint& virtDetPoint,
+                         const TYPE radius,
+                         Gaudi::XYZPoint& sphReflPoint ) const
+      {
+ 
+        // typedefs vectorised types
+        using Eigen3Vector = Eigen::Matrix< TYPE , 3 , 1 >;
+        using Vector = LHCb::Math::Eigen::XYZVector<TYPE>;
+        using Point  = LHCb::Math::Eigen::XYZPoint<TYPE>;
+        using Vec4x  = 
+          typename std::conditional<std::is_same<TYPE,float>::value,Vec4f,Vec4d>::type;
+
+        // vector from mirror centre of curvature to assumed emission point
+        const Vector evec( emissionPoint - CoC );
+        const TYPE e2 = evec.dot(evec);
+
+        // vector from mirror centre of curvature to virtual detection point
+        const Vector dvec( virtDetPoint - CoC );
+        const TYPE d2 = dvec.dot(dvec);
+
+        // various quantities needed to create quartic equation
+        // see LHCB/98-040 section 3, equation 3
+        const auto ed2 = e2 * d2;
+        const TYPE cosgamma2 = ( ed2 > 0 ? std::pow(evec.dot(dvec),2)/ed2 : 1.0 );
+        // vectorise 4 square roots into 1
+        const auto tmp_sqrt = sqrt( Vec4x( e2, d2,
+                                           cosgamma2 < 1.0 ? 1.0-cosgamma2 : 0.0,
+                                           cosgamma2 ) );
+        const auto e         = tmp_sqrt[0];
+        const auto d         = tmp_sqrt[1];
+        const auto singamma  = tmp_sqrt[2];
+        const auto cosgamma  = tmp_sqrt[3];
+
+        const auto dx        = d * cosgamma;
+        const auto dy        = d * singamma;
+        const auto r2        = radius * radius;
+        const auto dy2       = dy * dy;
+        const auto edx       = e + dx;
+
+        // Fill array for quartic equation
+        const auto a0      =     4.0 * ed2;
+        const auto inv_a0  =   ( a0 > 0 ? 1.0 / a0 : std::numeric_limits<TYPE>::max() );
+        const auto dyrad2  =     2.0 * dy * radius;
+        const auto a1      = - ( 2.0 * dyrad2 * e2 ) * inv_a0;
+        const auto a2      =   ( (dy2 * r2) + ( edx * edx * r2 ) - a0 ) * inv_a0;
+        const auto a3      =   ( dyrad2 * e * (e-dx) ) * inv_a0;
+        const auto a4      =   ( ( e2 - r2 ) * dy2 ) * inv_a0;
+
+        // use simplified RICH version of quartic solver
+        const auto sinbeta = solve_quartic_RICH<TYPE>( a1, // a
+                                                       a2, // b
+                                                       a3, // c
+                                                       a4  // d
+                                                       );
+
+        // (normalised) normal vector to reflection plane
+        auto n = evec.cross3(dvec);
+        n /= std::sqrt( n.dot(n) );
+
+        // construct rotation transformation
+        // Set vector magnitude to radius
+        // rotate vector and update reflection point
+        const Eigen::AngleAxis<TYPE> angleaxis( vdt::fast_asinf(sinbeta),
+                                                Eigen3Vector(n[0],n[1],n[2]) );
+        sphReflPoint = ( CoC + 
+                         Gaudi::XYZVector( angleaxis *
+                                           Eigen3Vector(evec[0],evec[1],evec[2]) *
+                                           ( radius / e ) ) );
+
+      }
+
     };
 
   }
 }
-
-#endif // RICHRECPHOTONTOOLS_QuarticSolver_H
diff --git a/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
new file mode 100755
index 0000000000000000000000000000000000000000..dfc5da3a9f3e2e8f388a6d1de300b0ad1acb6481
--- /dev/null
+++ b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
@@ -0,0 +1,324 @@
+
+//----------------------------------------------------------------------
+/** @file QuarticSolverNewton.h
+ *
+ *  @author Christina Quast, Rainer Schwemmer
+ *  @date   2017-02-03
+ */
+//----------------------------------------------------------------------
+
+#pragma once
+
+// Gaudi
+#include "GaudiKernel/Kernel.h"
+
+// VectorClass
+#include "VectorClass/vectorclass.h"
+
+// STL
+#include <math.h>
+#include <type_traits>
+#include <array>
+
+namespace Rich
+{
+  namespace Rec
+  {
+
+    //-----------------------------------------------------------------------------
+    /** @class QuarticSolverNewton
+     *
+     *  Utility class that implements the solving of the Quartic equation for the RICH
+     *  Based on original code by Chris Jones
+     *
+     *  @author Christina Quast, Rainer Schwemmer
+     *  @date   2017-02-03
+     */
+    //-----------------------------------------------------------------------------
+    class QuarticSolverNewton
+    {
+
+    public:
+
+      /** Solves the characteristic quartic equation for the RICH optical system.
+       *
+       *  See note LHCB/98-040 RICH section 3 for more details
+       *
+       *  @param emissionPoint Assumed photon emission point on track
+       *  @param CoC           Spherical mirror centre of curvature
+       *  @param virtDetPoint  Virtual detection point
+       *  @param radius        Spherical mirror radius of curvature
+       *  @param sphReflPoint  The reconstructed reflection pont on the spherical mirror
+       *
+       *  @return boolean indicating status of the quartic solution
+       *  @retval true  Calculation was successful. sphReflPoint is valid.
+       *  @retval false Calculation failed. sphReflPoint is not valid.
+       */
+      template< class TYPE,
+                std::size_t BISECTITS = 3, std::size_t NEWTONITS = 4,
+                class POINT = Gaudi::XYZPoint >
+      inline void solve( const POINT& emissionPoint,
+                         const POINT& CoC,
+                         const POINT& virtDetPoint,
+                         const TYPE radius,
+                         POINT& sphReflPoint ) const
+      {
+        using namespace std;
+
+        // vector from mirror centre of curvature to assumed emission point
+        const auto evec( emissionPoint - CoC );
+        const TYPE e2 = evec.Dot(evec);
+
+        // vector from mirror centre of curvature to virtual detection point
+        const auto dvec( virtDetPoint - CoC );
+        const TYPE d2 = dvec.Dot(dvec);
+
+        // various quantities needed to create quartic equation
+        // see LHCB/98-040 section 3, equation 3
+        const TYPE       ed2 = e2 * d2;
+        const TYPE cosgamma2 = ( ed2 > TYPE(0.0) ? pow(evec.Dot(dvec),2)/ed2 : TYPE(1.0) );
+        // vectorise 4 square roots into 1
+        using Vec4x = typename std::conditional<std::is_same<TYPE,float>::value,Vec4f,Vec4d>::type;
+        const auto tmp_sqrt = sqrt( Vec4x( e2, d2,
+                                           cosgamma2 < TYPE(1.0) ? TYPE(1.0)-cosgamma2 : TYPE(0.0),
+                                           cosgamma2 ) );
+        const TYPE e         = tmp_sqrt[0];
+        const TYPE d         = tmp_sqrt[1];
+        const TYPE singamma  = tmp_sqrt[2];
+        const TYPE cosgamma  = tmp_sqrt[3];
+        // scalar versions
+        //const TYPE e         = std::sqrt(e2);
+        //const TYPE d         = std::sqrt(d2);
+        //const TYPE singamma  = std::sqrt( cosgamma2 < TYPE(1.0) ? TYPE(1.0)-cosgamma2 : TYPE(0.0) );
+        //const TYPE cosgamma  = std::sqrt(cosgamma2);
+
+        const TYPE dx        = d * cosgamma;
+        const TYPE dy        = d * singamma;
+        const TYPE r2        = radius * radius;
+        const TYPE dy2       = dy * dy;
+        const TYPE edx       = e + dx;
+
+        // Fill array for quartic equation
+        const TYPE a0        =     TYPE(4.0) * ed2;
+        //Newton solver doesn't care about a0 being not 1.0. Remove costly division and several multiplies.
+        //This has some downsides though. The a-values are hovering around a numerical value of 10^15. single
+        //precision float max is 10^37. A single square and some multiplies will push it over the limit of what
+        //single precision float can handle. It's ok for the newton method, but Halley or higher order Housholder
+        //will fail without this normalization.
+        //const auto inv_a0  =   ( a0 > 0 ? 1.0 / a0 : std::numeric_limits<TYPE>::max() );
+        const TYPE dyrad2  =     TYPE(2.0) * dy * radius;
+        const TYPE a1      = - ( TYPE(2.0) * dyrad2 * e2 );// * inv_a0;
+        const TYPE a2      =   ( (dy2 * r2) + ( edx * edx * r2 ) - a0 );// * inv_a0;
+        const TYPE a3      =   ( dyrad2 * e * (e-dx) );// * inv_a0;
+        const TYPE a4      =   ( ( e2 - r2 ) * dy2 );// * inv_a0;
+
+        // use simplified RICH version of quartic solver
+        //        const auto sinbeta = solve_quartic_RICH<TYPE>( a1, // a
+        //                                                       a2, // b
+        //                                                       a3, // c
+        //                                                       a4  // d
+        //                                                       );
+
+        //Use optimized newton solver on quartic equation.
+        const auto sinbeta = solve_quartic_newton_RICH<TYPE,BISECTITS,NEWTONITS>( a0, a1, a2, a3, a4 );
+
+        //TODO: This method should be better but has problems still for some reasons
+        //const auto sinbeta = solve_quartic_housholder_RICH<TYPE, 3>(a1, a2, a3, a4);
+
+        // construct rotation transformation
+        // Set vector magnitude to radius
+        // rotate vector and update reflection point
+        //rotation matrix uses sin(beta) and cos(beta) to perform rotation
+        //even fast_asinf (which is only single precision and defeats the purpose
+        //of this class being templatable to double btw) is still too slow
+        //plus there is a cos and sin call inside AngleAxis ...
+        //We can do much better by just using the cos(beta) we already have to calculate
+        //sin(beta) and do our own rotation. On top of that we rotate non-normalized and save several
+        //Divisions by normalizing only once at the very end
+        //Again, care has to be taken since we are close to float_max here without immediate normalization.
+        //As far as we have tried with extreme values in the rich coordinate systems this is fine.
+        const TYPE nx  = (evec.y()*dvec.z()) - (evec.z()*dvec.y());
+        const TYPE ny  = (evec.z()*dvec.x()) - (evec.x()*dvec.z());
+        const TYPE nz  = (evec.x()*dvec.y()) - (evec.y()*dvec.x());
+
+        const TYPE nx2 = nx*nx; 
+        const TYPE ny2 = ny*ny; 
+        const TYPE nz2 = nz*nz; 
+
+        const TYPE norm      = nx2 + ny2 + nz2;
+        const TYPE norm_sqrt = sqrt(norm);
+
+        const TYPE        a = sinbeta * norm_sqrt;
+        const TYPE sinbeta2 = sinbeta * sinbeta;
+        const TYPE        b = ( sinbeta2 < TYPE(1.0) ? ( TYPE(1.0) - sqrt( TYPE(1.0) - sinbeta2 ) ) : TYPE(1.0) );
+        const TYPE    enorm = radius/(e*norm);
+
+        const TYPE bnxny = b*nx*ny;
+        const TYPE bnxnz = b*nx*nz;
+        const TYPE bnynz = b*ny*nz;
+
+        //Perform non-normalized rotation
+        const std::array<TYPE,9> M = 
+          { norm - b*(nz2+ny2), a*nz+bnxny,         -a*ny+bnxnz,
+            -a*nz+bnxny,       norm - b*(nx2+nz2),  a*nx+bnynz,
+            a*ny+bnxnz,       -a*nx+bnynz,         norm-b*(ny2+nx2) };
+ 
+        //re-normalize rotation and scale to radius in one step
+        const TYPE ex = enorm*(evec.x()*M[0]+evec.y()*M[3]+evec.z()*M[6]);
+        const TYPE ey = enorm*(evec.x()*M[1]+evec.y()*M[4]+evec.z()*M[7]);
+        const TYPE ez = enorm*(evec.x()*M[2]+evec.y()*M[5]+evec.z()*M[8]);
+
+        //sphReflPoint = ( CoC + Gaudi::XYZVector( ex, ey, ez ) );
+        sphReflPoint = { CoC.x() + ex, CoC.y() + ey, CoC.z() + ez };
+      }
+
+    private:
+
+      //A newton iteration solver for the Rich quartic equation
+      //Since the polynomial that is evaluated here is extremely constrained
+      //(root is in small interval, one root guaranteed), we can use a much more
+      //efficient approximation (which still has the same precision) instead of the
+      //full blown mathematically absolute correct method and still end up with
+      //usable results
+
+      template < class TYPE >
+      inline TYPE f4( const TYPE &a0, const TYPE &a1, const TYPE &a2, 
+                      const TYPE &a3, const TYPE &a4, TYPE &x ) const
+      {
+        //return (a0 * x*x*x*x + a1 * x*x*x + a2 * x*x + a3 * x + a4);
+        //A bit more FMA friendly
+        return ( ( ( ( ( ( ( a0 * x ) + a1 ) * x ) + a2 ) * x ) + a3 ) * x ) + a4;
+      }
+
+      template < class TYPE >
+      inline TYPE df4( const TYPE &a0, const TYPE &a1, const TYPE &a2, const TYPE &a3, TYPE &x ) const
+      {
+        //return (4.0f*a0 * x*x*x + 3.0f*a1 * x*x + 2.0f*a2 * x + a3);
+        return ( ( ( ( TYPE(4.0) * a0*x ) + TYPE(3.0) * a1 ) * x + ( TYPE(2.0) * a2 ) ) * x ) + a3;
+      }
+
+      /** Horner's method to evaluate the polynomial and its derivatives with as little math operations as
+       *  possible. We use a template here to allow the compiler to unroll the for loops and produce code
+       *  that is free from branches and optimized for the grade of polynomial and derivatives as necessary.
+       */
+      template < class TYPE, std::size_t ORDER = 4, std::size_t DIFFGRADE = 3 >
+      inline void evalPolyHorner( const TYPE (&a)[ORDER+1], TYPE (&res)[DIFFGRADE+1], const TYPE x ) const
+      {
+        for ( std::size_t i = 0; i <= DIFFGRADE; ++i )
+        {
+          res[i] = a[0];
+        }
+
+        for ( std::size_t j = 1; j <= ORDER; ++j )
+        {
+          res[0] = res[0] * x + a[j];
+          const int l = (ORDER - j) > DIFFGRADE ? DIFFGRADE : ORDER - j;
+          for ( int i = 1; i <= l ; ++i )
+          {
+            res[i] = res[i] * x + res[i-1];
+          }
+        }
+
+        //TODO: Check assembly to see if this is optimized away if DIFFGRADE is <= 2
+        TYPE l = 1.0;
+        for ( std::size_t i = 2; i <= DIFFGRADE; ++i )
+        {
+          l *= i;
+          res[i] = res[i] * l;
+        }
+
+      }
+
+      /** 3rd grade Housholder's method for iteratively finding the root of a function. Tests have shown that we seem to be
+       *  already too constrained in our polynomial and input data that we are not gaining anything over newton. In fact it seems
+       *  to make performance worse.
+       *  TODO: find out why performance of this is so bad. We might be missing something. It should converge much faster than newton.
+       */
+      template < class TYPE, std::size_t ITER = 3 >
+      inline TYPE householder(const TYPE (&a)[5], TYPE x0) const
+      {
+        TYPE res[4];
+        for ( std::size_t i = 0; i < ITER; ++i )
+        {
+          evalPolyHorner<TYPE, 4, 3>(a, res, x0);
+          x0 -= ((TYPE(6.0) * res[0]*res[1]*res[1] - TYPE(3.0) * res[0]*res[0]*res[2])/
+                 (TYPE(6.0) * res[1]*res[1]*res[1] - TYPE(6.0) * res[0]*res[1]*res[2] + res[0]*res[0]*res[3]));
+        }
+        return x0;
+      }
+
+      /** Use Housholder's method to solve the rich quartic equation. Important! If this function is used, the coefficients of the
+       *  polynomial have to be normalized. They are typically around 10^15 and floating point overflows will occur in here if these
+       *  large values are used.
+       */
+      template < class TYPE, std::size_t ITER = 3 >
+      inline TYPE solve_quartic_housholder_RICH( const TYPE& a0,
+                                                 const TYPE& a1,
+                                                 const TYPE& a2,
+                                                 const TYPE& a3) const
+      {
+        //Starting value for housholder iteration. Empirically values seem to be
+        //between 0 0and 0.4 so we chose the value in the middle as starting point
+        TYPE x0 = TYPE(0.2);
+        TYPE a[5] = { TYPE(1.0), a0, a1, a2, a3 };
+        return householder<TYPE, ITER>( a, x0 );
+      }
+
+      /** Newton-Rhapson method for calculating the root of the rich polynomial. It uses the bisection method in the beginning
+       *  to get close enough to the root to allow the second stage newton method to converge faster. After 4 iterations of newton
+       *  precision is as good as single precision floating point will get you. We have introduced a few tuning parameters like the
+       *  newton gain factor and a slightly skewed bisection division, which in this particular case help to speed things up.
+       *  TODO: Once we are happy with the number of newton and bisection iterations, this function should be templated to the number of
+       *  these iterations to allow loop unrolling and elimination of unnecessary branching
+       *  TODO: These tuning parameters have been found by low effort experimentation on random input data. A more detailed
+       *  study should be done with real data to find the best values
+       */
+      template < class TYPE, std::size_t BISECTITS = 3, std::size_t NEWTONITS = 4 >
+      inline TYPE solve_quartic_newton_RICH( const TYPE& a0,
+                                             const TYPE& a1,
+                                             const TYPE& a2,
+                                             const TYPE& a3,
+                                             const TYPE& a4 ) const
+      {
+        //Use N steps of bisection method to find starting point for newton
+        TYPE l(0.0);
+        TYPE u(0.5);
+        TYPE m(0.2); //We start a bit off center since the distribution of roots tends to be more to the left side
+        const TYPE a[5] = {a0, a1, a2, a3, a4};
+        TYPE res[2];
+        for ( std::size_t i = 0; i < BISECTITS; ++i )
+        {
+          const auto oppositeSign = std::signbit(f4(a0,a1,a2,a3,a4,m) * f4(a0,a1,a2,a3,a4,l));
+
+          l = oppositeSign ? l : m;
+          u = oppositeSign ? m : u;
+          //0.4 instead of 0.5 to speed up convergence. Most roots seem to be closer to 0 than to the extreme end
+          m = ( u + l ) * TYPE(0.4);
+        }
+
+        //Newton for the rest
+        // CRJ : why copy the value. original value is not needed later on...
+        //TYPE x = m;
+
+        //Most of the times we are approaching the root of the polynomial from one side
+        //and fall short by a certain fraction. This fraction seems to be around 1.04 of the
+        //quotient which is subtracted from x. By scaling it up, we take bigger steps towards
+        //the root and thus converge faster.
+        //TODO: study this factor more closely it's pure guesswork right now. We might get
+        //away with 3 iterations if we can find an exact value
+        const TYPE gain = 1.04;
+
+        for ( std::size_t i = 0; i < NEWTONITS; ++i )
+        {
+          evalPolyHorner<TYPE, 4, 1>( a, res, m );
+          //epsilon = f4(a0,a1,a2,a3,a4,x) / df4(a0,a1,a2,a3,x);
+          m -= gain * ( res[0] / res[1] );
+        }
+
+        return m;
+      }
+
+    };
+
+  }
+}
diff --git a/Tr/TrackFitEvent/doc/release.notes b/Tr/TrackFitEvent/doc/release.notes
index ebe90661e0e49d37af3d3b2cc72eaedb6d3a43a0..861836c30abee94a39ae7d826356a33652c15d5e 100755
--- a/Tr/TrackFitEvent/doc/release.notes
+++ b/Tr/TrackFitEvent/doc/release.notes
@@ -4,6 +4,9 @@
 ! Purpose     : Package for the tracking fitting-related classes
 !-----------------------------------------------------------------------------
 
+! 2016-12-14 - Jeroen van Tilburg
+ - Adapt to new FT geometry and numbering scheme
+
 !========================= TrackFitEvent v6r6 2016-01-28 =========================
 ! 2015-12-05 - Gerhard Raven
  - reduce mention of 'static' and 'std::auto_ptr' 
diff --git a/Tr/TrackFitEvent/src/FTMeasurement.cpp b/Tr/TrackFitEvent/src/FTMeasurement.cpp
index d18b80557421aa60fff41c250d7f996eaf8583f5..bb21204db05b325ee4078aae0a1adb8e43868d01 100644
--- a/Tr/TrackFitEvent/src/FTMeasurement.cpp
+++ b/Tr/TrackFitEvent/src/FTMeasurement.cpp
@@ -5,7 +5,7 @@
 
 // from Event
 #include "Event/FTCluster.h"
-#include "Kernel/LineTraj.h"
+
 // local
 #include "Event/FTMeasurement.h"
 
@@ -18,17 +18,13 @@ using namespace LHCb;
 //-----------------------------------------------------------------------------
 void FTMeasurement::init( const DeFTDetector& geom ) {
 
-  const DeFTFibreMat* ftMat = geom.findFibreMat( m_cluster.channelID() );
+  const DeFTMat* ftMat = geom.findMat( m_cluster.channelID() );
   m_detectorElement = ftMat;
-  DetectorSegment mySeg = ftMat->createDetSegment( m_cluster.channelID(), m_cluster.fraction() );
-  Gaudi::XYZPoint begPoint( mySeg.x( mySeg.yMin() ), mySeg.yMin(), mySeg.z( mySeg.yMin() ) );
-  Gaudi::XYZPoint endPoint( mySeg.x( mySeg.yMax() ), mySeg.yMax(), mySeg.z( mySeg.yMax() ) );
-
-  m_size = m_cluster.size();
-  m_errMeasure = 0.050 + .030 * m_size;
-  m_trajectory.reset( new LHCb::LineTraj( begPoint, endPoint ) );
-  m_z = ftMat->layerCenterZ();
-  m_measure = mySeg.x( 0. );
+  m_size = m_cluster.isLarge();
+  m_errMeasure = 0.080; // + .030 * m_size;
+  m_trajectory = ftMat->trajectory(m_cluster.channelID(), m_cluster.fraction());
+  m_z = ftMat->globalZ();
+  m_measure = 0. ; // JvT: I believe this is purely historical. Should remove it?
 
 }
 
diff --git a/Tr/TrackFitEvent/xml/FTMeasurement.xml b/Tr/TrackFitEvent/xml/FTMeasurement.xml
index 3a2bed3ab7646c0f5406febbf97f68cab9022c41..e1f3a31fb0be8f40130356b642ee553ed5a21a30 100644
--- a/Tr/TrackFitEvent/xml/FTMeasurement.xml
+++ b/Tr/TrackFitEvent/xml/FTMeasurement.xml
@@ -44,7 +44,7 @@
         </code>
       </method>
 
-      <method
+      <!--method
         type  = 'double'
         name  = 'charge'
         desc  = 'Retrieve the cluster charge'
@@ -52,7 +52,7 @@
         <code>
           return m_cluster.charge();
         </code>
-      </method>
+      </method-->
 
 
       <constructor
diff --git a/Tr/TrackMCTools/doc/release.notes b/Tr/TrackMCTools/doc/release.notes
index 4e14d1a8882064704099f58b1fa67c42d37fc2cd..d631e91644827b7335d29a0ce4f805710ef5db23 100755
--- a/Tr/TrackMCTools/doc/release.notes
+++ b/Tr/TrackMCTools/doc/release.notes
@@ -4,6 +4,10 @@
 ! Purpose     : package for tracking tools accessing MC information
 !-----------------------------------------------------------------------------
 
+! 2016-12-14 - Jeroen van Tilburg
+ - Adapt to new FT geometry and numbering scheme
+ - Changed location of Linker tables. 
+
 !========================= TrackMCTools v3r13 2016-03-18 =========================
 ! 2016-02-08 - Adam Davis
  - Add more debugging tools for PatLongLivedTracking
diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp b/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp
index 4bcfc474c3ce058df99a5b5753135aafe35bcfe9..c8f0f80999675197740e78a2a3dcd79a4b6b546e 100755
--- a/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp
+++ b/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp
@@ -10,7 +10,7 @@
 #include "Event/VPCluster.h"
 #include "Event/OTTime.h"
 #include "Event/MuonCoord.h"
-#include "Event/FTCluster.h"
+#include "Event/FTLiteCluster.h"
 
 
 DECLARE_TOOL_FACTORY( LHCbIDsToMCHits )
@@ -219,7 +219,7 @@ void LHCbIDsToMCHits::linkFT(const LHCbID& lhcbid, LinkMap& output) const
 {
   if (!m_configuredFT){
     m_configuredFT = true;
-    m_ftLinks = FTLinks( evtSvc(), msgSvc(),LHCb::FTClusterLocation::Default+m_endString);
+    m_ftLinks = FTLinks( evtSvc(), msgSvc(),LHCb::FTLiteClusterLocation::Default+m_endString);
     if (m_ftLinks.notFound()) {
       throw GaudiException("no FTLinker", "LHCbIDsToMCHits" ,
                            StatusCode::FAILURE);
diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCHits.h b/Tr/TrackMCTools/src/LHCbIDsToMCHits.h
index 75ce6ead745e648276c48b80f5514f2b2cd41540..506c1bc64da4eb345dad335f5ec2b8fcc24b9cc0 100755
--- a/Tr/TrackMCTools/src/LHCbIDsToMCHits.h
+++ b/Tr/TrackMCTools/src/LHCbIDsToMCHits.h
@@ -24,8 +24,7 @@ namespace LHCb{
   class VeloCluster;
   class MuonCoord;
   class VPCluster;
-  class UTCluster;
-  class FTCluster;
+  class UTCluster;  
 }
 
 class LHCbIDsToMCHits: public GaudiTool, virtual public ILHCbIDsToMCHits,
@@ -90,9 +89,7 @@ private:
   typedef LinkedTo<LHCb::MCHit,LHCb::VeloCluster> VeloLinks;
   typedef LinkedTo<LHCb::MCHit,LHCb::VPCluster> VPLinks;
   typedef LinkedTo<LHCb::MCHit,LHCb::MuonCoord> MuonLinks;
-  typedef LinkedTo<LHCb::MCHit,LHCb::FTCluster> FTLinks;
-
-
+  typedef LinkedTo<LHCb::MCHit> FTLinks;
 
   template<typename ID, typename LINKER>
   void linkToDetTruth(const ID& id, LINKER& aLinker, LinkMap& output ) const;
diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp
index 8a5a2eb761f10a28d7105a2ca006c7d6efddc730..d65fce775b6076e4f00814e1f33f9440b5f018b6 100755
--- a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp
+++ b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp
@@ -8,7 +8,7 @@
 #include "Event/STCluster.h"
 #include "Event/VeloCluster.h"
 #include "Event/VPCluster.h"
-#include "Event/FTCluster.h"
+#include "Event/FTLiteCluster.h"
 #include "Event/OTTime.h"
 #include "Event/MuonCoord.h"
 
@@ -237,7 +237,7 @@ StatusCode LHCbIDsToMCParticles::linkFT(const LHCbID& lhcbid, LinkMap& output) c
   if (!m_configuredFT){
     m_configuredFT = true;
     m_ftLinks = FTLinks(evtSvc(), msgSvc(),
-                            LHCb::FTClusterLocation::Default);
+                            LHCb::FTLiteClusterLocation::Default);
     if (m_ftLinks.notFound()) {
       return Error("no FTLinker",StatusCode::FAILURE,10);
     }
diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h
index 80718e92543e4875ff2915707640634fbd936dd7..8a6397df1fcfa7a54aa99369f2a5357e17ad7717 100755
--- a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h
+++ b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h
@@ -24,7 +24,6 @@ namespace LHCb{
   class VeloCluster;
   class MuonCoord;
   class VPCluster;
-  class FTCluster;
 }
 
 class LHCbIDsToMCParticles: public GaudiTool, virtual public ILHCbIDsToMCParticles,
@@ -90,7 +89,7 @@ private:
   typedef LinkedTo<LHCb::MCParticle,LHCb::MuonCoord> MuonLinks;
   // -- upgrade
   typedef LinkedTo<LHCb::MCParticle,LHCb::VPCluster> VPLinks;
-  typedef LinkedTo<LHCb::MCParticle,LHCb::FTCluster> FTLinks;
+  typedef LinkedTo<LHCb::MCParticle> FTLinks;
 
   template<typename ID, typename LINKER>
   void linkToDetTruth(const ID& id, LINKER& aLinker, LinkMap& output ) const;
diff --git a/Tr/TrackTools/doc/release.notes b/Tr/TrackTools/doc/release.notes
index 403a8e9a81dd3b4b817ccd368d497a8b9b0e74f1..944ec15f975577d41695b037faab8c235a7ee258 100755
--- a/Tr/TrackTools/doc/release.notes
+++ b/Tr/TrackTools/doc/release.notes
@@ -4,6 +4,10 @@
 ! Purpose     : Package for tracking tools not accessing MC information
 !-----------------------------------------------------------------------------
 
+! 2016-12-14 - Jeroen van Tilburg
+ - Adapt to new FT geometry and numbering scheme
+ - Changed location of Linker tables. 
+
 ! 2016-06-22 - David Hutchcroft
  - Added code to set the pT of the VELO only tracks with a sign based on the first strip 
    in the first cluster rather than the track key (which was random). Also an option to reverse the sign 
diff --git a/Tr/TrackTools/src/FTHitExpectation.cpp b/Tr/TrackTools/src/FTHitExpectation.cpp
index 18cd99f1e4fb23d443e571f0129b13500b31607a..671f3aef6f346d3929f04c070adf33ae2359b70c 100644
--- a/Tr/TrackTools/src/FTHitExpectation.cpp
+++ b/Tr/TrackTools/src/FTHitExpectation.cpp
@@ -5,18 +5,13 @@
 #include "Event/Track.h"
 #include "Event/State.h"
 
-//#include "Kernel/ISTChannelIDSelector.h"
-
 // Tsa
 #include "TsaKernel/Line.h"
 #include "TsaKernel/Parabola.h"
 #include "TsaKernel/Line3D.h"
-//#include "TsaKernel/IITExpectedHits.h"
 
 #include "LHCbMath/GeomFun.h"
 #include "FTDet/DeFTDetector.h"
-#include "FTDet/DeFTLayer.h"
-#include "FTDet/DeFTFibreMat.h"
 
 // local
 #include "FTHitExpectation.h"
@@ -36,30 +31,14 @@ DECLARE_TOOL_FACTORY( FTHitExpectation )
 //=============================================================================
 StatusCode FTHitExpectation::initialize()
 {
-
   StatusCode sc = THitExpectation::initialize();
-  if (sc.isFailure()){
-    return Error("Failed to initialize", sc);
-  }
+  if (sc.isFailure()) return Error("Failed to initialize", sc);
 
-  //m_expectedITHits = tool<Tf::Tsa::IITExpectedHits>("Tf::Tsa::ITExpectedHits");
-
-  // need to know the z of T stations....
+  // Get all layers (need to know the z-plane of each layer)
   m_ftDet = getDet<DeFTDetector>(DeFTDetectorLocation::Default);
-  if ( m_ftDet -> version() < 20 ){
-    return Error("FTDetector version not implemented (different than v2 or v5)",StatusCode::FAILURE);
-  }
-
-  m_layers    = m_ftDet -> layers();
-  m_fibremats = m_ftDet -> fibremats();
+  if ( m_ftDet -> version() < 61 )
+    return Error("This version requires FTDet v6.1 or higher", StatusCode::FAILURE);
 
-  /*
-  // (selector
-  if (m_selectorName != "ALL"){
-    info() << "Selector name " << m_selectorName << endmsg;
-    m_selector  = tool< ISTChannelIDSelector >( m_selectorType,m_selectorName );
-  }
-  */
   return StatusCode::SUCCESS;
 }
 
@@ -69,135 +48,100 @@ IHitExpectation::Info FTHitExpectation::expectation(const LHCb::Track& aTrack) c
   IHitExpectation::Info info;
   info.likelihood = 0.0;
   info.nFound = 0;
-  info.nExpected = 0;;
-
+  info.nExpected = 0;
 
   const auto& ids = aTrack.lhcbIDs();
   std::vector<LHCb::LHCbID> ftHits; ftHits.reserve(ids.size());
   std::copy_if( ids.begin(), ids.end(), std::back_inserter(ftHits),
                 [](const LHCb::LHCbID& id) { return id.isFT(); } );
 
-  for ( auto iterL = m_layers.begin(); iterL != m_layers.end() ; ++iterL )
-    {
-      LHCb::FTChannelID elemID( (unsigned int)std::distance( m_layers.begin(), iterL ), 0u, 0u, 0u, 0u );
-      double layerZ = ((*iterL)->geometry()->toGlobal( Gaudi::XYZPoint(0. ,0., 0.) )).z();
-
-      Tf::Tsa::Line     line    = yLine    ( aTrack, layerZ );
-      Tf::Tsa::Parabola aParab  = xParabola( aTrack, layerZ );
-      Tf::Tsa::Line tanLine     = aParab.tangent( layerZ );
-      Tf::Tsa::Line3D aLine3D   = Tf::Tsa::createLine3D( tanLine, line, layerZ );
-      std::vector<std::pair<LHCb::FTChannelID, double> > vectFTPairs;
-
-      collectHits( aLine3D, elemID, vectFTPairs );
-
-      int old = -1;
-      for ( const auto& p : vectFTPairs )
-	{
-	  if ( int(p.second) != old) {
-	    old = int(p.second);
-	    if ( p.first ) {
-	      ++(info.nExpected);
-	      auto itIter = std::find( ftHits.begin(), ftHits.end(), p.first);
-	      if (itIter != ftHits.end() ) ++info.nFound;
-	    }
-	  }
-	}  // pairs
-    }//layers
-  return info;
-}
+  // Loop over all layers
+  for( auto station : m_ftDet->stations() ) {
+    for( auto layer : station->layers() ) {
+
+      const DeFTMat* mat;
+      Gaudi::XYZPoint intersectionPoint;
+      bool expectHit = findIntersectingMat(layer, aTrack, mat, intersectionPoint);
+
+      if( expectHit ) {
+        ++(info.nExpected);
 
-void FTHitExpectation::collectHits( Tf::Tsa::Line3D aLine3D,
-				    LHCb::FTChannelID chan,
-				    std::vector<std::pair<LHCb::FTChannelID, double> >& vectFTPairs
-				    ) const{
-
-  for ( const auto& f : m_fibremats)
-    {
-      if( (unsigned int)(f -> layer()) != chan.layer() ) continue;
-      // propagate to z of sector
-      Gaudi::Plane3D entryplane, exitplane;
-      EntryExitPlanes( *f, entryplane, exitplane );
-
-      Gaudi::XYZPoint globalEntry = intersection( aLine3D, entryplane );
-      Gaudi::XYZPoint globalExit  = intersection( aLine3D, exitplane );
-
-      LHCb::MCHit hit;
-      hit.setEntry( globalEntry );
-      hit.setDisplacement ( globalExit-globalEntry  );
-      hit.setEnergy( 1.0 );
-      if( f -> isInside( hit.midPoint() ) ){
-        f -> calculateHits( hit, vectFTPairs ).ignore();
+        // Check of this layer is present in the hits
+        auto itIter = std::find_if( ftHits.begin(), ftHits.end(), [&mat] (const LHCb::LHCbID& id)
+            {return id.ftID().uniqueLayer() == mat->elementID().uniqueLayer(); } );
+        if (itIter != ftHits.end() ) ++info.nFound;
       }
     }
-}
-
+  }
 
-void FTHitExpectation::EntryExitPlanes( DeFTFibreMat& mat,  Gaudi::Plane3D& entry,  Gaudi::Plane3D& exit ) const
-{
-  double xmin = mat.fibreMatMinX();
-  double ymin = mat.fibreMatMinY();
-  double xmax = mat.fibreMatMaxX();
-  double ymax = mat.fibreMatMaxY();
-
-  double xLower = -0.5*(xmax-xmin);
-  double xUpper =  0.5*(xmax-xmin);
-  double yLower = -0.5*(ymax-ymin);
-  double yUpper =  0.5*(ymax-ymin);
-
-  const Gaudi::XYZPoint g1 = mat.geometry() -> toGlobal( Gaudi::XYZPoint( xLower, yLower, 0. ) );
-  const Gaudi::XYZPoint g2 = mat.geometry() -> toGlobal( Gaudi::XYZPoint( xLower, yUpper, 0. ) );
-  // trajectory of middle
-  //const Gaudi::XYZPoint g3 = mat.geometry()->toGlobal(xLower, 0., 0.);
-  const Gaudi::XYZPoint g4 = mat.geometry()->toGlobal(Gaudi::XYZPoint( xUpper, 0., 0.) );
-  // plane
-  Gaudi::Plane3D plane(g1,g2,g4);
-
-  double thickness = 1.2;
-  entry = Gaudi::Plane3D(plane.Normal(), mat.geometry()->toGlobal( Gaudi::XYZPoint(0.,0.,-0.5*thickness)));
-  exit  = Gaudi::Plane3D(plane.Normal(), mat.geometry()->toGlobal( Gaudi::XYZPoint(0.,0., 0.5*thickness)));
+  return info;
 }
 
 
-void FTHitExpectation::collect(const LHCb::Track& aTrack ,std::vector<LHCb::LHCbID>& ids) const{
+bool FTHitExpectation::findIntersectingMat(const DeFTLayer* layer, const LHCb::Track& aTrack,
+    const DeFTMat*& mat, Gaudi::XYZPoint& intersectionPoint) const {
+
+  // make a Line3D from the track
+  double layerZ = layer->globalZ();
+  Tf::Tsa::Line     line    = yLine    ( aTrack, layerZ );
+  Tf::Tsa::Parabola aParab  = xParabola( aTrack, layerZ );
+  Tf::Tsa::Line tanLine     = aParab.tangent( layerZ );
+  Tf::Tsa::Line3D aLine3D   = Tf::Tsa::createLine3D( tanLine, line, layerZ );
 
- for ( auto iterL = m_layers.begin(); iterL != m_layers.end() ; ++iterL )
-    {
-      LHCb::FTChannelID elemID( (unsigned int)std::distance( m_layers.begin(), iterL ), 0u, 0u, 0u, 0u );
-      double layerZ = ((*iterL)->geometry()->toGlobal( Gaudi::XYZPoint(0. ,0., 0.) )).z();
+  // find intersection point of track and plane of layer
+  const Gaudi::Plane3D layerPlane = layer->plane();
+  intersectionPoint = intersection(aLine3D, layerPlane);
 
-      Tf::Tsa::Line     line    = yLine    ( aTrack, layerZ );
-      Tf::Tsa::Parabola aParab  = xParabola( aTrack, layerZ );
-      Tf::Tsa::Line tanLine     = aParab.tangent( layerZ );
-      Tf::Tsa::Line3D aLine3D   = Tf::Tsa::createLine3D( tanLine, line, layerZ );
-      std::vector<std::pair<LHCb::FTChannelID, double> > vectFTPairs;
+  // find the module that is hit
+  const DeFTModule* module = layer->findModule( intersectionPoint );
+  if( module == nullptr ) return false;
 
-      collectHits( aLine3D, elemID, vectFTPairs );
+  // find intersection point of track and plane of module
+  // (to account for misalignments between module and layer)
+  tanLine     = aParab.tangent( intersectionPoint.z() );
+  aLine3D   = Tf::Tsa::createLine3D( tanLine, line, intersectionPoint.z() );
+  const Gaudi::Plane3D modulePlane = module->plane();
+  intersectionPoint = intersection(aLine3D, modulePlane);
 
+  // Find the corresponding mat
+  mat = module->findMat(intersectionPoint);
 
-      int old = -1;
-      ids.reserve(vectFTPairs.size());
-      for ( std::vector<std::pair<LHCb::FTChannelID, double> >::iterator iter = vectFTPairs.begin();
-	    iter != vectFTPairs.end();
-	    ++iter )
-	{
-	  if ( int(iter->second) != old ) {
-	    old = int(iter->second);
-	    if ( iter->first) ids.push_back( LHCb::LHCbID(iter->first) );
-	  }
-	}  // pairs
-    }//layers
+  // check if intersection point is inside the fibres of the module
+  return (mat == nullptr ) ? false : mat->isInside( intersectionPoint );
 }
 
+void FTHitExpectation::collect(const LHCb::Track& aTrack,
+                               std::vector<LHCb::LHCbID>& ids) const{
+  // Loop over all layers
+  for( auto station : m_ftDet->stations() ) {
+    for( auto layer : station->layers() ) {
+
+      const DeFTMat* mat;
+      Gaudi::XYZPoint intersectionPoint;
+      bool expectHit = findIntersectingMat(layer, aTrack, mat, intersectionPoint);
+
+      if( expectHit ) {
+        // Find the channel that is closest
+        Gaudi::XYZPoint localP = mat->geometry()->toLocal( intersectionPoint );
+        double frac;
+        LHCb::FTChannelID ftChan = mat->calculateChannelAndFrac(localP.x(), frac);
+
+        // Add the channels
+        // JvT: Without the fraction the bare FTChannelID is pretty useless...
+        if( std::abs(frac) <= 0.5 ) ids.push_back( LHCb::LHCbID(ftChan) );
+      }
+    }
+  }
+}
 
 
 unsigned int FTHitExpectation::nExpected(const LHCb::Track& aTrack) const{
-
   return expectation(aTrack).nExpected;
 }
 
+
 Gaudi::XYZPoint FTHitExpectation::intersection(const Tf::Tsa::Line3D& line,
 					       const Gaudi::Plane3D& aPlane) const{
-  // make a plane
   Gaudi::XYZPoint inter;
   double mu = 0;
   Gaudi::Math::intersection( line, aPlane, inter, mu );
diff --git a/Tr/TrackTools/src/FTHitExpectation.h b/Tr/TrackTools/src/FTHitExpectation.h
index 010f9fcad1e23f050d12961d2ff3c972f2b09bf9..2fe8ba3de7034c709f48f3f4344a0fb0e7d42521 100644
--- a/Tr/TrackTools/src/FTHitExpectation.h
+++ b/Tr/TrackTools/src/FTHitExpectation.h
@@ -47,7 +47,6 @@ public:
   */
   unsigned int nExpected ( const LHCb::Track& aTrack ) const override;
 
-
   /** Returns number of hits expected
   *
   *  @param aTrack Reference to the Track to test
@@ -67,27 +66,13 @@ public:
 
 private:
 
+  bool findIntersectingMat(const DeFTLayer* layer, const LHCb::Track& aTrack,
+      const DeFTMat*& mat, Gaudi::XYZPoint& intersectionPoint) const ;
 
-  //bool select(const LHCb::FTChannelID& chan) const;
-
-  void EntryExitPlanes( DeFTFibreMat& mat,  Gaudi::Plane3D& entry,  Gaudi::Plane3D& exit ) const;
-  void collectHits( Tf::Tsa::Line3D aLine3D,
-		    LHCb::FTChannelID chan,
-		    std::vector<std::pair<LHCb::FTChannelID, double> >& vectFTPairs
-		    ) const;
   Gaudi::XYZPoint intersection(const Tf::Tsa::Line3D& line, const Gaudi::Plane3D& aPlane) const;
 
-
   DeFTDetector* m_ftDet = nullptr;
-  std::vector<DeFTLayer*>    m_layers;
-  std::vector<DeFTFibreMat*> m_fibremats;
-
+ 
 };
 
-/*
-inline bool ITHitExpectation::select(const LHCb::STChannelID& chan) const{
-  return !m_selector || m_selector->select(chan);
-}
-*/
-
 #endif
diff --git a/Tr/TrackTools/src/FTMeasurementProvider.cpp b/Tr/TrackTools/src/FTMeasurementProvider.cpp
index f3beec1cd9c7a05f9d6d9ab07be503a8575afc67..12978959fa04794a10a0cb2de9ec338c246a6a28 100644
--- a/Tr/TrackTools/src/FTMeasurementProvider.cpp
+++ b/Tr/TrackTools/src/FTMeasurementProvider.cpp
@@ -31,8 +31,11 @@ public:
   using base_class::base_class;
 
   StatusCode initialize() override;
-  LHCb::Measurement* measurement( const LHCb::LHCbID& id, bool localY=false ) const  override;
-  LHCb::Measurement* measurement( const LHCb::LHCbID& id, const LHCb::ZTrajectory& refvector, bool localY=false) const override;
+  LHCb::Measurement* measurement( const LHCb::LHCbID& id,
+                                  bool localY=false ) const override;
+  LHCb::Measurement* measurement( const LHCb::LHCbID& id,
+                                  const LHCb::ZTrajectory& refvector,
+                                  bool localY=false) const override;
   inline LHCb::FTMeasurement* ftmeasurement( const LHCb::LHCbID& id ) const  ;
 
   void handle ( const Incident& incident ) override;
@@ -91,13 +94,10 @@ const FastClusterContainer<LHCb::FTLiteCluster,int>* FTMeasurementProvider::clus
 {
   /// If there is a new event, get the clusters and sort them according to their channel ID
   if( UNLIKELY(!m_clusters) ){
-    m_clusters = getIfExists<FastClusterContainer<LHCb::FTLiteCluster,int> >( LHCb::FTLiteClusterLocation::Default );
-    if(m_clusters){ //TODO: sort before putting clusters on the TES!
-      std::sort(m_clusters->begin(), m_clusters->end(),
-                [](const LHCb::FTLiteCluster& lhs, const LHCb::FTLiteCluster& rhs){ return lhs.channelID() < rhs.channelID(); });
-    }else{
-      error() << "Could not find FTLiteClusters at: " <<  LHCb::FTLiteClusterLocation::Default << endmsg;
-    }
+    m_clusters = getIfExists<FastClusterContainer<LHCb::FTLiteCluster,int> >
+                 ( LHCb::FTLiteClusterLocation::Default );
+    if(!m_clusters) error() << "Could not find FTLiteClusters at: "
+                            << LHCb::FTLiteClusterLocation::Default << endmsg;
   }
   return m_clusters ;
 }
@@ -111,8 +111,8 @@ LHCb::FTMeasurement* FTMeasurementProvider::ftmeasurement( const LHCb::LHCbID& i
     error() << "Not an FT measurement" << endmsg ;
     return nullptr;
   }
-  /// The clusters are sorted, so we can use a binary search (lower bound search) to find the element
-  /// corresponding to the channel ID
+  /// The clusters are sorted, so we can use a binary search (lower bound search)
+  /// to find the element corresponding to the channel ID
   const auto& c = *clusters();
   auto itH = std::lower_bound( c.begin(),  c.end(), id.ftID(),
                       [](const LHCb::FTLiteCluster clus, const LHCb::FTChannelID id){
@@ -124,7 +124,8 @@ LHCb::FTMeasurement* FTMeasurementProvider::ftmeasurement( const LHCb::LHCbID& i
 //-----------------------------------------------------------------------------
 /// Return the measurement
 //-----------------------------------------------------------------------------
-LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id, bool /*localY*/ ) const {
+LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id,
+                                                       bool /*localY*/ ) const {
   return ftmeasurement(id) ;
 }
 
@@ -133,7 +134,8 @@ LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id, b
 //-----------------------------------------------------------------------------
 
 LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id,
-                                                       const LHCb::ZTrajectory& /* reftraj */, bool /*localY*/ ) const {
+                                                       const LHCb::ZTrajectory& /* reftraj */,
+                                                       bool /*localY*/ ) const {
   // default implementation
   return ftmeasurement(id) ;
 }