diff --git a/Trigger/TrigAccel/TrigAccelEvent/TrigAccelEvent/TrigITkAccelEDM.h b/Trigger/TrigAccel/TrigAccelEvent/TrigAccelEvent/TrigITkAccelEDM.h
index 2663e9f24238b1d5f50b95782b695b7ca1421464..4ffa354aa35359237a704845361bb35db74c6189 100644
--- a/Trigger/TrigAccel/TrigAccelEvent/TrigAccelEvent/TrigITkAccelEDM.h
+++ b/Trigger/TrigAccel/TrigAccelEvent/TrigAccelEvent/TrigITkAccelEDM.h
@@ -37,6 +37,7 @@ namespace ITk {
     int m_nLayers;
     int m_nModules;
     SILICON_LAYER m_layers[MAX_SILICON_LAYERS];
+    int m_middleSpacePointLayers[MAX_SILICON_LAYERS];
     int m_hashArray[MAX_NUMBER_PIX_MODULES+MAX_NUMBER_SCT_MODULES];
     float m_minRZ[MAX_NUMBER_PIX_MODULES+MAX_NUMBER_SCT_MODULES];
     float m_maxRZ[MAX_NUMBER_PIX_MODULES+MAX_NUMBER_SCT_MODULES];
@@ -53,6 +54,7 @@ namespace ITk {
     int m_nSpacepoints;
     int m_nPhiSlices;
     int m_nLayers;
+    int m_nMiddleLayers;
     int m_index[MAX_NUMBER_SPACEPOINTS];
     int m_type[MAX_NUMBER_SPACEPOINTS];
     float m_x[MAX_NUMBER_SPACEPOINTS];
@@ -62,6 +64,7 @@ namespace ITk {
     float m_phi[MAX_NUMBER_SPACEPOINTS];
     float m_covR[MAX_NUMBER_SPACEPOINTS];
     float m_covZ[MAX_NUMBER_SPACEPOINTS];
+    float m_clusterWidth[MAX_NUMBER_SPACEPOINTS];
     SPACEPOINT_LAYER_RANGE m_phiSlices[MAX_PHI_SLICES];
   } SPACEPOINT_STORAGE;
 
diff --git a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/TrigInDetAccelerationService/ITrigInDetAccelerationSvc.h b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/TrigInDetAccelerationService/ITrigInDetAccelerationSvc.h
index be940bba0cdf3c1bb5890cd26ee1d393dd6d98ed..55745294a1dd150908fabc27391d315e9e1ebc15 100644
--- a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/TrigInDetAccelerationService/ITrigInDetAccelerationSvc.h
+++ b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/TrigInDetAccelerationService/ITrigInDetAccelerationSvc.h
@@ -31,6 +31,7 @@ class ITrigInDetAccelerationSvc : virtual public IService {
   //helper
 
   virtual const std::vector<short>& getLayerInformation(int) const = 0;
+  virtual size_t getMiddleLayersSize() const = 0;
 
 };
 
diff --git a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.cxx b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.cxx
index 9e8dea50e3774fd3fd0884154f90a48dabf52544..a25bab74eedb83d04169435bcd5a14795694428b 100644
--- a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.cxx
+++ b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.cxx
@@ -36,6 +36,7 @@ TrigInDetAccelerationSvc::TrigInDetAccelerationSvc( const std::string& name, ISv
   declareProperty( "NumberOfDCs", m_nDCs = 8 );
   declareProperty( "ModuleName", m_moduleName = "libTrigInDetCUDA.so");
   declareProperty( "useITkGeometry", m_useITkGeometry = false);
+  declareProperty( "MiddleSpacePointLayers", m_middleSpacePointLayers = std::vector<int>(), "Global IDs of layers that can contain middle spacepoints of track seeds" );
 } 
 
 
@@ -211,6 +212,7 @@ bool TrigInDetAccelerationSvc::exportITkGeometryInformation(const std::map<std::
   pArray->m_nModules=0;
   
   int layerIdx=0;
+  int middleLayerIdx=0;
 
   for(std::map<std::tuple<short,short, int, int>,std::vector<PhiEtaHash> >::const_iterator it = hashMap.begin();it!=hashMap.end();++it, layerIdx++) {
     
@@ -223,6 +225,12 @@ bool TrigInDetAccelerationSvc::exportITkGeometryInformation(const std::map<std::
     pArray->m_layers[layerIdx].m_nElements = 0;
     pArray->m_layers[layerIdx].m_subdet = globalLayerId;
     pArray->m_layers[layerIdx].m_type = barrel_ec;
+
+    // Fill in a table of layer ids that can contain the middle SPs
+    if (std::find(m_middleSpacePointLayers.begin(), m_middleSpacePointLayers.end(), globalLayerId) != m_middleSpacePointLayers.end()){
+      pArray->m_middleSpacePointLayers[middleLayerIdx] = layerIdx;
+      ++middleLayerIdx;
+    }
     
     std::vector<std::vector<PhiEtaHash>::const_iterator> vStops;
     vStops.push_back((*it).second.begin());
diff --git a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.h b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.h
index 8780aef36874637aba3170bda56aa2e265f84850..4b765af181a4749d9f153347178b5514b71cef72 100644
--- a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.h
+++ b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationService/src/TrigInDetAccelerationSvc.h
@@ -44,6 +44,7 @@ class TrigInDetAccelerationSvc : public extends<AthService, ITrigInDetAccelerati
   
   virtual TrigAccel::Work* createWork(unsigned int, std::shared_ptr<TrigAccel::OffloadBuffer>) const override;
   virtual const std::vector<short>& getLayerInformation(int) const override;
+  virtual size_t getMiddleLayersSize() const override {return m_middleSpacePointLayers.size();};
 
  private:   
 
@@ -70,6 +71,7 @@ class TrigInDetAccelerationSvc : public extends<AthService, ITrigInDetAccelerati
   int m_nDCs;
   std::string m_moduleName;
   bool m_useITkGeometry;
+  std::vector<int> m_middleSpacePointLayers;
   void* m_libHandle; //for OffloadFactory
   TrigAccel::WorkFactory* m_pWF;
   TrigAccel::Module* m_module;
diff --git a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx
index 338371c3e33ec7ee671debeb9dbb1ebf3f87f2e7..a700381a82fc6769e604eb681f318503c45eeaa8 100644
--- a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx
+++ b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx
@@ -6,6 +6,7 @@
 #include "AthenaBaseComps/AthMsgStreamMacros.h" 
 #include "AthenaBaseComps/AthCheckMacros.h"
 #include "TrigAccelEvent/TrigITkAccelEDM.h"
+#include "InDetPrepRawData/PixelCluster.h"
 
 TrigITkAccelerationTool::TrigITkAccelerationTool(const std::string& t, 
 						     const std::string& n,
@@ -83,7 +84,9 @@ size_t TrigITkAccelerationTool::exportSeedMakingJob(const TrigCombinatorialSetti
   sfs.m_maxTripletBufferLength = tcs.m_maxTripletBufferLength;
   sfs.m_isFullScan = 1;
 
-  if(!(roi->isFullscan() || roi->composite() )){
+  // roi->isFullscan() check cannot be currently used for standalone ITk ROI configuration
+  bool isRoiFullscan = ((roi->etaPlus() >= 4) && (roi->etaMinus() <= -4));
+  if(!(isRoiFullscan || roi->composite() )){
     //roi suitable for gpu
     //composite rois are not supported at this point
     sfs.m_isFullScan = 0;
@@ -146,14 +149,21 @@ size_t TrigITkAccelerationTool::exportSeedMakingJob(const TrigCombinatorialSetti
   sps.m_nSpacepoints = nSP;
   sps.m_nPhiSlices = nSlices;
   sps.m_nLayers = nLayers;
+  sps.m_nMiddleLayers = m_accelSvc->getMiddleLayersSize();
 
   int spIdx=0;
   for(int slice = 0;slice<nSlices;slice++) {
     for(int layer = 0;layer<nLayers;layer++) {
       int layerStart = spIdx;
+      bool isBarrel = (layerTypes[layer] == 0);
       std::vector<std::pair<int, const TrigSiSpacePointBase*> >& v = phiLArray[layer + slice*nLayers];
       for(std::vector<std::pair<int, const TrigSiSpacePointBase*> >::iterator it = v.begin();it!=v.end();++it) {
         const TrigSiSpacePointBase* sp  = (*it).second;
+        const Trk::SpacePoint* osp = sp->offlineSpacePoint();
+        const InDet::PixelCluster* pCL = dynamic_cast<const InDet::PixelCluster*>(osp->clusterList().first);
+        float clusterWidth = pCL->width().widthPhiRZ().y();
+        if (!isBarrel && clusterWidth > 0.2) continue;
+
         sps.m_index[spIdx] = (*it).first;
         sps.m_type[spIdx] = 1; // always Pixel
         sps.m_x[spIdx] = sp->x();
@@ -163,6 +173,7 @@ size_t TrigITkAccelerationTool::exportSeedMakingJob(const TrigCombinatorialSetti
         sps.m_phi[spIdx] = sp->phi();
         sps.m_covR[spIdx] = sp->dr()*sp->dr();
         sps.m_covZ[spIdx] = sp->dz()*sp->dz();
+        sps.m_clusterWidth[spIdx] = clusterWidth;
         spIdx++;
       }
       int layerEnd = spIdx;
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletCountingKernelCuda_ITk.cuh b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletCountingKernelCuda_ITk.cuh
index b0f5b6714316c8a5377a082d27ca5f54e0c83f15..9528be7910bf8e24ffa546838603b7cfa9dd6565 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletCountingKernelCuda_ITk.cuh
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletCountingKernelCuda_ITk.cuh
@@ -30,7 +30,10 @@ __global__ static void doubletCountingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 	const float maxOuterRadius = 550.0;
 	 
 	const int sliceIdx = blockIdx.x;
-	const int layerIdx = blockIdx.y;
+	const int layerIdx = dDetModel->m_middleSpacePointLayers[blockIdx.y];
+
+	const TrigAccel::ITk::SILICON_LAYER& layerGeo =  dDetModel->m_layers[layerIdx];
+	bool isBarrel = (layerGeo.m_type == 0);
 
 	if(threadIdx.x == 0 && threadIdx.y == 0) {
 		const TrigAccel::ITk::SPACEPOINT_LAYER_RANGE& slr = dSpacepoints->m_phiSlices[sliceIdx];
@@ -64,7 +67,14 @@ __global__ static void doubletCountingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 		float zm = dSpacepoints->m_z[spmIdx];
 		float rm = dSpacepoints->m_r[spmIdx];
 
-		if (!canBeMiddleSpacePoint(rm)) continue;
+		float minTau = 0;
+		float maxTau = 100;
+
+		if (isBarrel) {
+			float clusterWidth = dSpacepoints->m_clusterWidth[spmIdx];
+			minTau = 6.7*(clusterWidth - 0.2);
+			maxTau = 1.6 + 0.15/(clusterWidth + 0.2) + 6.1*(clusterWidth - 0.2);
+		}
 
 		if(threadIdx.y ==0) {
 			nInner[threadIdx.x] = 0;
@@ -92,16 +102,16 @@ __global__ static void doubletCountingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 				if(next_spEnd == next_spBegin) continue;//no spacepoints in this layer
 
 				// Check if SPs from the nextLayer (currently processed layer) and the current SP would exceed the maximum doublet length
-				const TrigAccel::ITk::SILICON_LAYER& layerGeo =  dDetModel->m_layers[nextLayerIdx];
-				bool isBarrel = (layerGeo.m_type == 0);
-				float refCoord = layerGeo.m_refCoord; // Position of the layer in r (barrel) or z (end cap)
-				if(isBarrel && std::abs(refCoord-rm)>maxDoubletLength) continue;
+				const TrigAccel::ITk::SILICON_LAYER& layerGeoSp =  dDetModel->m_layers[nextLayerIdx];
+				bool isBarrelSp = (layerGeoSp.m_type == 0);
+				float refCoord = layerGeoSp.m_refCoord; // Position of the layer in r (barrel) or z (end cap)
+				if(isBarrelSp && std::abs(refCoord-rm)>maxDoubletLength) continue;
 
 				// Calculate the nextLayer boundaries by projecting the current spacepoint within the beamspot limits
 				float minCoord = 10000.0;
 				float maxCoord =-10000.0;
 
-				if(isBarrel) {
+				if(isBarrelSp) {
 					minCoord = zMinus + refCoord*(zm-zMinus)/rm;
 					maxCoord = zPlus + refCoord*(zm-zPlus)/rm;
 				}
@@ -114,7 +124,7 @@ __global__ static void doubletCountingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 					float tmp = maxCoord;maxCoord = minCoord;minCoord = tmp;
 				}
 
-				if(layerGeo.m_maxBound<minCoord || layerGeo.m_minBound>maxCoord) continue;
+				if(layerGeoSp.m_maxBound<minCoord || layerGeoSp.m_minBound>maxCoord) continue;
 
 				//3. get a tile of inner/outer spacepoints
 
@@ -124,7 +134,7 @@ __global__ static void doubletCountingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 					float rsp = dSpacepoints->m_r[spIdx];
 
 					// Check if SP is within the doublet boundaries
-					float spCoord = (isBarrel) ? zsp : rsp;
+					float spCoord = (isBarrelSp) ? zsp : rsp;
 					if(spCoord<minCoord || spCoord>maxCoord) continue;
 
 					float dr = rsp - rm;
@@ -132,12 +142,18 @@ __global__ static void doubletCountingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 
 					// Cut on doublet length
 					float dL = std::sqrt( dr*dr + dz*dz);
-					float maxDL = getMaxDeltaLEta(getEta(dr, dz, dL));
+					float maxDL = GPUTrackSeedingItkHelpers::getMaxDeltaLEta(GPUTrackSeedingItkHelpers::getEta(dr, dz, dL));
 					if(std::abs(dL)>maxDL || std::abs(dL)<minDoubletLength) continue;
 								
-					// Cut on tau
+					// Cut on tau = cot(theta)
 					float tau = dz/dr;
-					if(std::abs(tau)>maxCtg) continue;
+					float ftau = std::abs(tau);
+					if(ftau>maxCtg) continue;
+
+					// Cut on pixel width
+					if (isBarrelSp) {
+						if(ftau < minTau || ftau > maxTau) continue;
+					}
 
 					// Cut on Z
 					float z0 = zsp - rsp*tau;
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh
index 12a09a10aee41906a1c1f9942d10e3d7cf136049..40cd9b2804998d10d8f8d30af22b389fd02e2d35 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh
@@ -11,6 +11,8 @@
  */
 
 
+namespace GPUTrackSeedingItkHelpers {
+
 /**
  * @brief Calculate eta for a doublet
  * @param dr radius difference between doublet's space points
@@ -30,10 +32,14 @@ __device__ static float getMaxDeltaLEta (float eta) {
   else return eta*eta*eta*eta*1.7582417 + eta*eta*-129.67033 + 3324.61538;
 }
 
+/**
+ * @brief Calculate |dt| cut for a given significance
+ * @param s std::sqrt(dt2/(covdt+dCov))
+ */
+__device__ static float getSignificanceCut (float s) {
+  return -0.0598 * std::log(s - 0.0664) - 0.0121;
+}
 
-__device__ static float canBeMiddleSpacePoint (float r) {
-  if(std::abs(r - 100.0) > 20.0) return false;
-  return true;
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMakingKernelCuda_ITk.cuh b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMakingKernelCuda_ITk.cuh
index 5e5d97060913fad43817f7805acf7f88231d1471..1f07c696424cbddaafba09c87a705c59f428c649 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMakingKernelCuda_ITk.cuh
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMakingKernelCuda_ITk.cuh
@@ -33,7 +33,10 @@ __global__ static void doubletMakingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SETTI
 	const float maxOuterRadius = 550.0;
 	
 	const int sliceIdx = blockIdx.x;
-	const int layerIdx = blockIdx.y;
+	const int layerIdx = dDetModel->m_middleSpacePointLayers[blockIdx.y];
+
+	const TrigAccel::ITk::SILICON_LAYER& layerGeo = dDetModel->m_layers[layerIdx];
+	bool isBarrel = (layerGeo.m_type == 0);
 
 	if(threadIdx.x == 0 && threadIdx.y == 0) {
 		const TrigAccel::ITk::SPACEPOINT_LAYER_RANGE& slr = dSpacepoints->m_phiSlices[sliceIdx];
@@ -67,7 +70,7 @@ __global__ static void doubletMakingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SETTI
 	for(int spmIdx=threadIdx.x+spBegin;spmIdx<spEnd;spmIdx+=blockDim.x) {
 
 		if(threadIdx.y ==0) {
-			 hasDoublets = d_Info->m_good[spmIdx] == 1;
+			hasDoublets = d_Info->m_good[spmIdx] == 1;
 		}
 		__syncthreads();
 
@@ -91,7 +94,14 @@ __global__ static void doubletMakingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SETTI
 		float zm = dSpacepoints->m_z[spmIdx];
 		float rm = dSpacepoints->m_r[spmIdx];
 
-		if (!canBeMiddleSpacePoint(rm)) continue;
+		float minTau = 0;
+		float maxTau = 100;
+
+		if (isBarrel) {
+			float clusterWidth = dSpacepoints->m_clusterWidth[spmIdx];
+			minTau = 6.7*(clusterWidth - 0.2);
+			maxTau = 1.6 + 0.15/(clusterWidth + 0.2) + 6.1*(clusterWidth - 0.2);
+		}
 
 		//2. loop over other phi-bins / layers
 
@@ -110,18 +120,18 @@ __global__ static void doubletMakingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SETTI
 
 				if(next_spEnd == next_spBegin) continue;//no spacepoints in this layer
 
-				const TrigAccel::ITk::SILICON_LAYER& layerGeo =  dDetModel->m_layers[nextLayerIdx];
-				bool isBarrel = (layerGeo.m_type == 0);
+				const TrigAccel::ITk::SILICON_LAYER& layerGeoSp =  dDetModel->m_layers[nextLayerIdx];
+				bool isBarrelSp = (layerGeoSp.m_type == 0);
 
-				float refCoord = layerGeo.m_refCoord;
-				if(isBarrel && std::abs(refCoord-rm)>maxDoubletLength) continue;
+				float refCoord = layerGeoSp.m_refCoord;
+				if(isBarrelSp && std::abs(refCoord-rm)>maxDoubletLength) continue;
 
 				//boundaries for nextLayer
 
 				float minCoord = 10000.0;
 				float maxCoord =-10000.0;
 
-				if(isBarrel) {
+				if(isBarrelSp) {
 					minCoord = zMinus + refCoord*(zm-zMinus)/rm;
 					maxCoord = zPlus + refCoord*(zm-zPlus)/rm;
 				}
@@ -134,7 +144,7 @@ __global__ static void doubletMakingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SETTI
 					float tmp = maxCoord;maxCoord = minCoord;minCoord = tmp;
 				}
 
-				if(layerGeo.m_maxBound<minCoord || layerGeo.m_minBound>maxCoord) continue;
+				if(layerGeoSp.m_maxBound<minCoord || layerGeoSp.m_minBound>maxCoord) continue;
 
 				//3. get a tile of inner/outer spacepoints
 
@@ -143,7 +153,7 @@ __global__ static void doubletMakingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SETTI
 					float zsp = dSpacepoints->m_z[spIdx];
 					float rsp = dSpacepoints->m_r[spIdx];
 
-					float spCoord = (isBarrel) ? zsp : rsp;
+					float spCoord = (isBarrelSp) ? zsp : rsp;
 
 					if(spCoord<minCoord || spCoord>maxCoord) continue;
 
@@ -152,12 +162,18 @@ __global__ static void doubletMakingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SETTI
 
 					// Cut on doublet length
 					float dL = std::sqrt( dr*dr + dz*dz);
-					float maxDL = getMaxDeltaLEta(getEta(dr, dz, dL));
+					float maxDL = GPUTrackSeedingItkHelpers::getMaxDeltaLEta(GPUTrackSeedingItkHelpers::getEta(dr, dz, dL));
 					if(std::abs(dL)>maxDL || std::abs(dL)<minDoubletLength) continue;
 								
 					// Cut on tau
 					float tau = dz/dr;
-					if(std::abs(tau)>maxCtg) continue;
+					float ftau = std::abs(tau);
+					if(ftau>maxCtg) continue;
+
+					// Cut on pixel width
+					if (isBarrelSp) {
+						if(ftau < minTau || ftau > maxTau) continue;
+					}
 
 					// Cut on Z
 					float z0 = zsp - rsp*tau;
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh
index b234c0e6ba6795998b62a1a31746f8b0e34e297c..a2f17efe924973bcdbf829d36ae12cc90db103a7 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh
@@ -38,21 +38,15 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 	__shared__ float cosA;
 	__shared__ float sinA;
 
-	//  __shared__ float R2inv_array[MAX_NUMBER_DOUBLETS];
-	__shared__ float Rinv_array[MAX_NUMBER_DOUBLETS_ITk];
-	__shared__ float t_array[MAX_NUMBER_DOUBLETS_ITk];
+	__shared__ float Rinv_array[MAX_NUMBER_DOUBLETS_ITk]; // inverse radius (xy plane)
+	__shared__ float tau_array[MAX_NUMBER_DOUBLETS_ITk]; // tau = cot(theta) = dz/dr
 	__shared__ int spIdx_array[MAX_NUMBER_DOUBLETS_ITk];
 	__shared__ float u_array[MAX_NUMBER_DOUBLETS_ITk];
 	__shared__ float v_array[MAX_NUMBER_DOUBLETS_ITk];
-
-	//  __shared__ float covZ_array[MAX_NUMBER_DOUBLETS];
-	// __shared__ float covR_array[MAX_NUMBER_DOUBLETS];
-
-	__shared__ float tCov_array[MAX_NUMBER_DOUBLETS_ITk];
-
+	__shared__ float tauCov_array[MAX_NUMBER_DOUBLETS_ITk]; // covariance of tau
 
 	__shared__ int PairIdx_array[MAX_TRIPLETS_ITk];
-	__shared__ float Q_array[MAX_TRIPLETS_ITk];
+	__shared__ float Q_array[MAX_TRIPLETS_ITk]; // Quality score for a triplet Q=d0*d0
 	__shared__ int sortedIdx[MAX_TRIPLETS_ITk];
 
 
@@ -61,7 +55,7 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 	__shared__ int nPairs;
 	__shared__ int nTriplets;
 
-	const double dtCut = 0.25;
+	const float dtCut = 0.3; // Cut on cot(theta) difference between two doublets
 	const float radLen = 0.036;
 	const float dp = 13.6/dSettings->m_tripletPtMin;
 	const float CovMS = dp*dp*radLen;
@@ -73,8 +67,8 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 	const float maxD0 = dSettings->m_tripletD0Max;
 
 	const float phiPlus = dSettings->m_phiPlus;
-  const float phiMinus = dSettings->m_phiMinus;
-  const bool isFullscan = (dSettings->m_isFullScan == 1);
+	const float phiMinus = dSettings->m_phiMinus;
+	const bool isFullscan = (dSettings->m_isFullScan == 1);
 
 
 	for(int itemIdx = blockIdx.x;itemIdx<maxItem;itemIdx += gridDim.x) {
@@ -109,9 +103,6 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 		}
 		__syncthreads();
 
-		if (!canBeMiddleSpacePoint(rm)) continue;    
-
-		
 		for(int innerIdx = threadIdx.x; innerIdx<nInner;innerIdx+=blockDim.x) {
 		 
 			int k = atomicAdd(&iDoublet,1);  
@@ -128,9 +119,9 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 
 				float R2inv = 1.0/(dx_inn*dx_inn+dy_inn*dy_inn); 
 				Rinv_array[k] = sqrt(R2inv);	
-				t_array[k] = Rinv_array[k]*dz_inn;
+				tau_array[k] = Rinv_array[k]*dz_inn;
 
-				tCov_array[k] = R2inv*(covZ + dSpacepoints->m_covZ[spiIdx] + t_array[k]*t_array[k]*(covR + dSpacepoints->m_covR[spiIdx]));
+				tauCov_array[k] = R2inv*(covZ + dSpacepoints->m_covZ[spiIdx] + tau_array[k]*tau_array[k]*(covR + dSpacepoints->m_covR[spiIdx]));
 				
 				float xn_inn = dx_inn*cosA + dy_inn*sinA; 
 				float yn_inn =-dx_inn*sinA + dy_inn*cosA;	
@@ -165,9 +156,9 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 
 				float R2inv = 1.0/(dx_out*dx_out+dy_out*dy_out);
 				Rinv_array[k] = sqrt(R2inv);
-				t_array[k] = Rinv_array[k]*dz_out;
+				tau_array[k] = Rinv_array[k]*dz_out;
 
-				tCov_array[k] = R2inv*(covZ + dSpacepoints->m_covZ[spoIdx] + t_array[k]*t_array[k]*(covR + dSpacepoints->m_covR[spoIdx]));
+				tauCov_array[k] = R2inv*(covZ + dSpacepoints->m_covZ[spoIdx] + tau_array[k]*tau_array[k]*(covR + dSpacepoints->m_covR[spoIdx]));
 
 				float xn_out = dx_out*cosA + dy_out*sinA; 
 				float yn_out =-dx_out*sinA + dy_out*cosA;	
@@ -192,23 +183,24 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 		//retrieve shared data for doublets doublet_i and doublet_j and apply cut(s)	  
 		
 		//0. dt matching
-		float t_inn = t_array[doublet_i];
-		float t_out = t_array[doublet_j];
-		float dt = t_inn - t_out;	
-		if(std::abs(dt)>dtCut) continue;
+		float tau_inn = tau_array[doublet_i];
+		float tau_out = tau_array[doublet_j];
+		float dt = tau_inn - tau_out;	
+		if (std::abs(dt)>dtCut) continue;
 		
 		//1. rz matching
 
-		float t_inn2 = t_inn*t_inn;
-		float tCov_inn = tCov_array[doublet_i];
-		float tCov_out = tCov_array[doublet_j];
+		float tau_inn2 = tau_inn*tau_inn;
+		float tauCov_inn = tauCov_array[doublet_i];
+		float tauCov_out = tauCov_array[doublet_j];
 
-		double dCov = CovMS*(1+t_inn2);
+		double dCov = CovMS*(1+tau_inn2);
 		
-		float covdt = tCov_inn + tCov_out; 
-		covdt += 2*Rinv_array[doublet_i]*Rinv_array[doublet_j]*(t_inn*t_out*covR + covZ); 
+		float covdt = tauCov_inn + tauCov_out; 
+		covdt += 2*Rinv_array[doublet_i]*Rinv_array[doublet_j]*(tau_inn*tau_out*covR + covZ); 
 		float dt2 = dt*dt*(1/9.0);
 		if(dt2 > covdt+dCov) continue;//i.e. 3-sigma cut 
+		if (GPUTrackSeedingItkHelpers::getSignificanceCut(std::sqrt(dt2/(covdt+dCov))) < std::abs(dt)) continue;
 
 		//2. pT estimate 
 
@@ -253,7 +245,6 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
 		}
 
 		//Calculate Quality    
-
 		float Q = d0*d0;
 
 		int l = atomicAdd(&nTriplets, 1);
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingWorkCuda_ITk.cu b/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingWorkCuda_ITk.cu
index 03645ebb4ce3ebef922de88a2740d005433e045d..f4882eb1428ce77d0ba8f2488fadfd9cb8d3aab0 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingWorkCuda_ITk.cu
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingWorkCuda_ITk.cu
@@ -212,23 +212,15 @@ bool SeedMakingWorkCudaManagedITk::run() {
   TrigAccel::ITk::OUTPUT_SEED_STORAGE* ps = reinterpret_cast<TrigAccel::ITk::OUTPUT_SEED_STORAGE*>(p.m_outputseeds);
   
   cudaSetDevice(id);
-
   checkError();
   
-
-
   cudaMemPrefetchAsync(p.m_settings, sizeof(TrigAccel::ITk::SEED_FINDER_SETTINGS), id, p.m_stream);
-
   checkError();
 
-
-
   cudaMemPrefetchAsync(p.m_spacepoints, sizeof(TrigAccel::ITk::SPACEPOINT_STORAGE), id, p.m_stream);
-
   checkError();
 
   cudaStreamSynchronize(p.m_stream);
-    
   TrigAccel::ITk::SEED_FINDER_SETTINGS* dSettings  = reinterpret_cast<TrigAccel::ITk::SEED_FINDER_SETTINGS *>(p.m_settings);
   TrigAccel::ITk::SPACEPOINT_STORAGE* dSpacepoints = reinterpret_cast<TrigAccel::ITk::SPACEPOINT_STORAGE *>(p.m_spacepoints);
   TrigAccel::ITk::DETECTOR_MODEL* dDetModel        = reinterpret_cast<TrigAccel::ITk::DETECTOR_MODEL*>(p.d_detmodel);
@@ -238,21 +230,20 @@ bool SeedMakingWorkCudaManagedITk::run() {
   DOUBLET_STORAGE_ITk* dStorage                   = reinterpret_cast<DOUBLET_STORAGE_ITk*>(p.d_doubletstorage);
 
   cudaMemset(p.m_outputseeds,0,10*sizeof(int));
-
   checkError();
 
   cudaMemset(p.d_doubletstorage,0,3*sizeof(int));
-
   checkError();
   
   const TrigAccel::ITk::SPACEPOINT_STORAGE* pSPS = reinterpret_cast<const TrigAccel::ITk::SPACEPOINT_STORAGE *>(p.m_spacepoints);
   int nSlices = pSPS->m_nPhiSlices;
+  int nMiddleLayers = pSPS->m_nMiddleLayers; // The first two kernels are launched for middle spacepoints
   int nLayers = pSPS->m_nLayers;
   
   int nMiddleSp = NUM_MIDDLE_THREADS_ITk;//determines size of the doublet/triplet buffers
   int nOtherSp = OUTER_THREADS_MULTIPLIER_ITk*p.m_gpuParams.m_nNUM_SMX_CORES/NUM_MIDDLE_THREADS_ITk;//the size of the spacepoint buffer
 
-  dim3 gridDimensions(nSlices, nLayers);
+  dim3 gridDimensions(nSlices, nMiddleLayers);
   dim3 blockDimensions(nMiddleSp, nOtherSp);
 
   cudaMemset(p.d_doubletinfo,0,sizeof(DOUBLET_INFO_ITk));
@@ -264,14 +255,12 @@ bool SeedMakingWorkCudaManagedITk::run() {
   checkError();
 
   doubletCountingKernel_ITk<<<gridDimensions, blockDimensions, 0, p.m_stream>>>(dSettings, dSpacepoints, dDetModel, dInfo, nLayers, nSlices);
-
   cudaStreamSynchronize(p.m_stream);
 
   checkError();
 
   doubletMakingKernel_ITk<<<gridDimensions, blockDimensions, 0, p.m_stream>>>(dSettings, dSpacepoints, dDetModel, dOutput, 
     dInfo, dStorage, nLayers, nSlices);
-
   cudaStreamSynchronize(p.m_stream);
 
   checkError();
@@ -283,7 +272,6 @@ bool SeedMakingWorkCudaManagedITk::run() {
   
   doubletMatchingKernel_ITk<<<p.m_gpuParams.m_nNUM_TRIPLET_BLOCKS, NUM_TRIPLET_THREADS_ITk, 0, p.m_stream>>>(dSettings, dSpacepoints, dDetModel, dInfo, 
     dStorage,  dOutput, nStats[0]);
-
   cudaStreamSynchronize(p.m_stream);
 
   checkError();
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/TrigITkModuleCuda.cu b/Trigger/TrigAccel/TrigInDetCUDA/src/TrigITkModuleCuda.cu
index 4e8acc6561026f2ddd793f9a256c63d786f15169..d43f0ce2a75b9829c1476021a80e9ef17096b57c 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/TrigITkModuleCuda.cu
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/TrigITkModuleCuda.cu
@@ -175,7 +175,7 @@ TrigAccel::Work* TrigITkModuleCuda::createWork(int workType, std::shared_ptr<Tri
     int deviceId = 0;//always using device 0 for the time being
 
     cudaSetDevice(deviceId);
-    cudaMalloc(&d_detmodel, sizeof(TrigAccel::DETECTOR_MODEL));
+    cudaMalloc(&d_detmodel, sizeof(TrigAccel::ITk::DETECTOR_MODEL));
     checkError();
 
     m_d_detmodel_ptrs[deviceId] = d_detmodel;
diff --git a/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py b/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py
index a2c0afcb3f6f84b92594ace0524c3142b1ee7f95..b626d1ecdf0127d113717e44adcff72c98435e69 100644
--- a/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py
+++ b/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py
@@ -21,6 +21,12 @@ def ITkFastTrackFinderStandaloneCfg(flags, SiSPSeededTrackCollectionKey = None):
     
     inDetAccelSvc = CompFactory.TrigInDetAccelerationSvc("TrigInDetAccelerationSvc")
     inDetAccelSvc.useITkGeometry = True # Allows to read and export the ITk geometry
+    inDetAccelSvc.MiddleSpacePointLayers = [81000, 82000, 
+        91005, 90014, 92000, 92001, 92002, 92003, 92004, 92005, 92006, 92007, 92008, 92009, 92010,  
+        92011, 92012, 92013, 92014, 92015, 92016, 92017, 92018, 92019, 92020, 92021, 92022, 
+        71005, 70014, 72000, 72001, 72002, 72003, 72004, 72005, 72006, 72007, 72008, 72009, 72010, 
+        72011, 72012, 72013, 72014, 72015, 72016, 72017, 72018, 72019, 72020, 72021, 72022
+    ]
     acc.addService(inDetAccelSvc)
 
     acc.addPublicTool(CompFactory.TrigITkAccelerationTool(name = "TrigITkAccelerationTool_FTF"))