[DO NOT MERGE] remove duplicate pixels

DataFormats/SiPixelDigi/interface/SiPixelDigiConstants.h

@@ -55,7 +55,10 @@ namespace sipixelconstants {
     inline constexpr uint32_t getRow(uint32_t ww) { return ((ww >> ROW_shift) & ROW_mask); }
     inline constexpr uint32_t getDCol(uint32_t ww) { return ((ww >> DCOL_shift) & DCOL_mask); }
     inline constexpr uint32_t getPxId(uint32_t ww) { return ((ww >> PXID_shift) & PXID_mask); }
-  }  // namespace functions
+    inline constexpr uint32_t removeADC(uint32_t ww) {
+      return (ww >> ADC_bits);
+    }  // ADC_shift ==0: let's keep it simple
+  }    // namespace functions
 }  // namespace sipixelconstants
 
 #endif  // DataFormats_SiPixelDigi_interface_SiPixelDigiConstants

RecoLocalTracker/SiPixelClusterizer/plugins/SiPixelRawToClusterGPUKernel.cu

@@ -347,6 +347,7 @@ namespace pixelgpudetails {
     //if (threadIdx.x==0) printf("Event: %u blockIdx.x: %u start: %u end: %u\n", eventno, blockIdx.x, begin, end);
 
     int32_t first = threadIdx.x + blockIdx.x * blockDim.x;
+    int32_t lastWord = wordCounter - 1;
     for (int32_t iloop = first, nend = wordCounter; iloop < nend; iloop += blockDim.x * gridDim.x) {
       auto gIndex = iloop;
       xx[gIndex] = 0;
@@ -440,6 +441,16 @@ namespace pixelgpudetails {
         }
       }
 
+      // remove duplicate pixels (for the time being keep second to reproduce current CPU behaviour)
+      auto noADC = sipixelconstants::removeADC(ww);
+      // auto noADCm1 =  sipixelconstants::removeADC(gIndex==0 ? 0 : word[gIndex-1]);
+      auto noADCp1 = sipixelconstants::removeADC(gIndex == lastWord ? 0 : word[gIndex + 1]);
+      if (noADC == noADCp1) {
+        // auto globalPix = frameConversion(barrel, side, layer, rocIdInDetUnit, localPix);
+        // printf("dup pix at %d %d %d\n",detId.moduleId,globalPix.row,globalPix.col);
+        continue;
+      }
+
       pixelgpudetails::Pixel globalPix = frameConversion(barrel, side, layer, rocIdInDetUnit, localPix);
       xx[gIndex] = globalPix.row;  // origin shifting by 1 0-159
       yy[gIndex] = globalPix.col;  // origin shifting by 1 0-415
@@ -656,7 +667,8 @@ namespace pixelgpudetails {
       std::cout << "CUDA findClus kernel launch with " << blocks << " blocks of " << threadsPerBlock << " threads\n";
 #endif
 
-      findClus<false><<<blocks, threadsPerBlock, 0, stream>>>(digis_d.view().moduleInd(),
+      findClus<false><<<blocks, threadsPerBlock, 0, stream>>>(digis_d.view().rawIdArr(),
+                                                              digis_d.view().moduleInd(),
                                                               digis_d.view().xx(),
                                                               digis_d.view().yy(),
                                                               clusters_d.moduleStart(),
@@ -763,7 +775,8 @@ namespace pixelgpudetails {
     threadsPerBlock = 256;
     blocks = phase2PixelTopology::numberOfModules;
 
-    findClus<true><<<blocks, threadsPerBlock, 0, stream>>>(digis_d.view().moduleInd(),
+    findClus<true><<<blocks, threadsPerBlock, 0, stream>>>(digis_d.view().rawIdArr(),
+                                                           digis_d.view().moduleInd(),
                                                            digis_d.view().xx(),
                                                            digis_d.view().yy(),
                                                            clusters_d.moduleStart(),

RecoLocalTracker/SiPixelClusterizer/plugins/gpuClusterChargeCut.h

@@ -39,6 +39,8 @@ namespace gpuClustering {
     auto endModule = moduleStart[0];
     for (auto module = firstModule; module < endModule; module += gridDim.x) {
       auto firstPixel = moduleStart[1 + module];
+      while (id[firstPixel] == invalidModuleId)
+        ++firstPixel;  // could be duplicates!
       auto thisModuleId = id[firstPixel];
       assert(thisModuleId < nMaxModules);
       assert(thisModuleId == moduleId[module]);

RecoLocalTracker/SiPixelClusterizer/plugins/gpuClustering.h

@@ -39,7 +39,8 @@ namespace gpuClustering {
   }
 
   template <bool isPhase2>
-  __global__ void findClus(uint16_t const* __restrict__ id,           // module id of each pixel
+  __global__ void findClus(uint32_t* __restrict__ rawIdArr,
+                           uint16_t* __restrict__ id,                 // module id of each pixel
                            uint16_t const* __restrict__ x,            // local coordinates of each pixel
                            uint16_t const* __restrict__ y,            //
                            uint32_t const* __restrict__ moduleStart,  // index of the first pixel of each module
@@ -103,6 +104,14 @@ namespace gpuClustering {
           printf("too many pixels in module %d: %d > %d\n", thisModuleId, msize - firstPixel, maxPixInModule);
           msize = maxPixInModule + firstPixel;
         }
+        /*
+        /// dump
+       if (thisModuleId==1) {
+        printf ("dump for mod 1\n");
+        for (int k = int(firstPixel); k< msize; ++k) printf ("%d %d\n",x[k],y[k]);
+        printf ("end dump for mod 1\n");
+       }
+       */
       }
 
       __syncthreads();
@@ -114,6 +123,24 @@ namespace gpuClustering {
       __syncthreads();
 #endif
 
+      // remove duplicate
+      if (msize > 1)
+        for (int i = first; i < msize - 1; i += blockDim.x) {
+          if (id[i] == invalidModuleId)  // skip invalid pixels
+            continue;
+          for (int j = i + 1; j < msize; ++j) {
+            if (id[j] == invalidModuleId)  // skip invalid pixels
+              continue;
+            if (y[i] == y[j] && x[i] == x[j]) {
+              //  printf("found dup %d %d %d %d %d\n",i,j,id[i],x[i], y[i]);
+              id[i] = invalidModuleId;
+              rawIdArr[i] = 0;
+              break;
+            }
+          }
+        }
+      __syncthreads();
+
       // fill histo
       for (int i = first; i < msize; i += blockDim.x) {
         if (id[i] == invalidModuleId)  // skip invalid pixels
@@ -204,6 +231,7 @@ namespace gpuClustering {
             continue;
           auto l = nnn[k]++;
           assert(l < maxNeighbours);
+          // if (l>=5) printf("too many Neighbours! %d %d %d\n",thisModuleId, x[i], y[i]);
           nn[k][l] = *p;
         }
       }

RecoLocalTracker/SiPixelClusterizer/test/gpuClustering_t.h

@@ -31,13 +31,15 @@ int main(void) {
   constexpr SiPixelClusterThresholds clusterThresholds(kSiPixelClusterThresholdsDefaultPhase1);
 
   // these in reality are already on GPU
+  auto h_raw = std::make_unique<uint32_t[]>(numElements);
   auto h_id = std::make_unique<uint16_t[]>(numElements);
   auto h_x = std::make_unique<uint16_t[]>(numElements);
   auto h_y = std::make_unique<uint16_t[]>(numElements);
   auto h_adc = std::make_unique<uint16_t[]>(numElements);
   auto h_clus = std::make_unique<int[]>(numElements);
 
 #ifdef __CUDACC__
+  auto d_raw = cms::cuda::make_device_unique<uint32_t[]>(numElements, nullptr);
   auto d_id = cms::cuda::make_device_unique<uint16_t[]>(numElements, nullptr);
   auto d_x = cms::cuda::make_device_unique<uint16_t[]>(numElements, nullptr);
   auto d_y = cms::cuda::make_device_unique<uint16_t[]>(numElements, nullptr);
@@ -265,6 +267,7 @@ int main(void) {
 
     cms::cuda::launch(findClus<false>,
                       {blocksPerGrid, threadsPerBlock},
+                      d_raw.get(),
                       d_id.get(),
                       d_x.get(),
                       d_y.get(),
@@ -305,8 +308,15 @@ int main(void) {
     h_moduleStart[0] = nModules;
     countModules<false>(h_id.get(), h_moduleStart.get(), h_clus.get(), n);
     memset(h_clusInModule.get(), 0, maxNumModules * sizeof(uint32_t));
-    findClus<false>(
-        h_id.get(), h_x.get(), h_y.get(), h_moduleStart.get(), h_clusInModule.get(), h_moduleId.get(), h_clus.get(), n);
+    findClus<false>(h_raw.get(),
+                    h_id.get(),
+                    h_x.get(),
+                    h_y.get(),
+                    h_moduleStart.get(),
+                    h_clusInModule.get(),
+                    h_moduleId.get(),
+                    h_clus.get(),
+                    n);
 
     nModules = h_moduleStart[0];
     auto nclus = h_clusInModule.get();