Follow up to the Alpaka integration in CMSSW

CondFormats/SiPixelObjects/BuildFile.xml

@@ -1,20 +1,24 @@
 <use name="alpaka"/>
-<use name="DataFormats/DetId"/>
-<use name="DataFormats/SiPixelDetId"/>
 <use name="CalibFormats/SiPixelObjects"/>
-<use name="FWCore/Utilities"/>
-<use name="DataFormats/FEDRawData"/>
-<use name="CondFormats/Serialization"/>
+<use name="CondCore/DBOutputService"/>
 <use name="CondFormats/DataRecord"/>
-<use name="CondFormats/SiStripObjects"/>
-<use name="Geometry/TrackerGeometryBuilder"/>
 <use name="CondFormats/External"/>
+<use name="CondFormats/Serialization"/>
+<use name="CondFormats/SiStripObjects"/>
+<use name="DataFormats/DetId"/>
+<use name="DataFormats/FEDRawData"/>
 <use name="DataFormats/GeometryVector"/>
+<use name="DataFormats/Portable"/>
+<use name="DataFormats/SiPixelDetId"/>
+<use name="DataFormats/SoATemplate"/>
+<use name="DataFormats/TrackerCommon"/>
 <use name="FWCore/MessageLogger"/>
+<use name="FWCore/Utilities"/>
 <use name="Geometry/CommonDetUnit"/>
 <use name="Geometry/CommonTopologies"/>
-<use name="HeterogeneousCore/AlpakaInterface"/>
+<use name="Geometry/TrackerGeometryBuilder"/>
 <use name="HeterogeneousCore/AlpakaCore"/>
+<use name="HeterogeneousCore/AlpakaInterface"/>
 <flags ALPAKA_BACKENDS="1"/>
 <export>
   <lib name="1"/>

CondFormats/SiPixelObjects/interface/SiPixelMappingHost.h

@@ -1,10 +1,9 @@
-#ifndef CondFormats_SiPixelObjects_SiPixelMappingHost_h
-#define CondFormats_SiPixelObjects_SiPixelMappingHost_h
+#ifndef CondFormats_SiPixelObjects_interface_SiPixelMappingHost_h
+#define CondFormats_SiPixelObjects_interface_SiPixelMappingHost_h
 
-#include <alpaka/alpaka.hpp>
-#include "DataFormats/Portable/interface/PortableHostCollection.h"
 #include "CondFormats/SiPixelObjects/interface/SiPixelMappingLayout.h"
+#include "DataFormats/Portable/interface/PortableHostCollection.h"
 
 using SiPixelMappingHost = PortableHostCollection<SiPixelMappingSoA>;
 
-#endif  // CondFormats_SiPixelObjects_SiPixelMappingHost_h
+#endif  // CondFormats_SiPixelObjects_interface_SiPixelMappingHost_h

CondFormats/SiPixelObjects/interface/alpaka/SiPixelMappingDevice.h

@@ -2,11 +2,10 @@
 #define CondFormats_SiPixelObjects_interface_alpaka_SiPixelMappingDevice_h
 
 #include <cstdint>
-#include <alpaka/alpaka.hpp>
+
 #include "DataFormats/Portable/interface/alpaka/PortableCollection.h"
 #include "CondFormats/SiPixelObjects/interface/SiPixelMappingLayout.h"
-#include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESProducer.h"
-#include "DataFormats/Portable/interface/PortableHostCollection.h"
+#include "HeterogeneousCore/AlpakaInterface/interface/config.h"
 
 namespace ALPAKA_ACCELERATOR_NAMESPACE {
 

CondFormats/SiPixelObjects/test/BuildFile.xml

@@ -1,6 +1,10 @@
 <use name="CondCore/DBOutputService"/>
 <use name="CondFormats/DataRecord"/>
 <use name="CondFormats/SiPixelObjects"/>
+<use name="FWCore/Framework"/>
+<use name="FWCore/ParameterSet"/>
+<use name="FWCore/ServiceRegistry"/>
+<use name="Geometry/Records"/>
 <use name="Geometry/TrackerGeometryBuilder"/>
 <library name="SiPixelObjectsTest" file="*.cc">
   <flags EDM_PLUGIN="1"/>

Geometry/CommonTopologies/interface/SimplePixelTopology.h

@@ -214,26 +214,26 @@ namespace phase2PixelTopology {
 
   HOST_DEVICE_CONSTANT uint8_t layerPairs[2 * nPairs] = {
 
-      0,  1,  0,  4,  0,  16,  //BPIX1 (3)
-      1,  2,  1,  4,  1,  16,  //BPIX2 (6)
-      2,  3,  2,  4,  2,  16,  //BPIX3 & Forward (9)
+      0,  1,  0,  4,  0,  16,  // BPIX1 (3)
+      1,  2,  1,  4,  1,  16,  // BPIX2 (6)
+      2,  3,  2,  4,  2,  16,  // BPIX3 & Forward (9)
 
-      4,  5,  5,  6,  6,  7,  7,  8,  8,  9,  9,  10, 10, 11,  //POS (16)
-      16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23,  //NEG (23)
+      4,  5,  5,  6,  6,  7,  7,  8,  8,  9,  9,  10, 10, 11,  // POS (16)
+      16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23,  // NEG (23)
 
       0,  2,  0,  5,  0,  17, 0,  6,  0,  18,  // BPIX1 Jump (28)
       1,  3,  1,  5,  1,  17, 1,  6,  1,  18,  // BPIX2 Jump (33)
 
-      11, 12, 12, 13, 13, 14, 14, 15,  //Late POS (37)
-      23, 24, 24, 25, 25, 26, 26, 27,  //Late NEG (41)
+      11, 12, 12, 13, 13, 14, 14, 15,  // Late POS (37)
+      23, 24, 24, 25, 25, 26, 26, 27,  // Late NEG (41)
 
-      4,  6,  5,  7,  6,  8,  7,  9,  8,  10, 9,  11, 10, 12,  //POS Jump (48)
-      16, 18, 17, 19, 18, 20, 19, 21, 20, 22, 21, 23, 22, 24,  //NEG Jump (55)
+      4,  6,  5,  7,  6,  8,  7,  9,  8,  10, 9,  11, 10, 12,  // POS Jump (48)
+      16, 18, 17, 19, 18, 20, 19, 21, 20, 22, 21, 23, 22, 24,  // NEG Jump (55)
   };
   HOST_DEVICE_CONSTANT uint32_t layerStart[numberOfLayers + 1] = {0,
                                                                   108,
                                                                   324,
-                                                                  504,  //Barrel
+                                                                  504,  // Barrel
                                                                   756,
                                                                   864,
                                                                   972,
@@ -245,7 +245,7 @@ namespace phase2PixelTopology {
                                                                   1620,
                                                                   1796,
                                                                   1972,
-                                                                  2148,  //Fp
+                                                                  2148,  // Fp
                                                                   2324,
                                                                   2432,
                                                                   2540,
@@ -257,7 +257,7 @@ namespace phase2PixelTopology {
                                                                   3188,
                                                                   3364,
                                                                   3540,
-                                                                  3716,  //Np
+                                                                  3716,  // Np
                                                                   numberOfModules};
 
   HOST_DEVICE_CONSTANT int16_t phicuts[nPairs]{
@@ -332,7 +332,7 @@ namespace pixelTopology {
     static constexpr uint32_t maxCellsPerHit = 256;
     static constexpr uint32_t avgTracksPerHit = 10;
     static constexpr uint32_t maxNumberOfTuples = 256 * 1024;
-    //this is well above thanks to maxNumberOfTuples
+    // this is well above thanks to maxNumberOfTuples
     static constexpr uint32_t maxHitsForContainers = avgHitsPerTrack * maxNumberOfTuples;
     static constexpr uint32_t maxNumberOfDoublets = 5 * 512 * 1024;
     static constexpr uint32_t maxNumOfActiveDoublets = maxNumberOfDoublets / 8;
@@ -381,6 +381,7 @@ namespace pixelTopology {
 
     static constexpr uint16_t numberOfModules = 3892;
 
+    // 1024 bins, 10 bits
     static constexpr uint16_t clusterBinning = 1024;
     static constexpr uint16_t clusterBits = 10;
 
@@ -391,7 +392,7 @@ namespace pixelTopology {
     static constexpr uint16_t firstEndcapPos = 4;
     static constexpr uint16_t firstEndcapNeg = 16;
 
-    static constexpr int16_t xOffset = -1e4;  //not used actually, to suppress static analyzer warnings
+    static constexpr int16_t xOffset = -1e4;  // not used actually, to suppress static analyzer warnings
 
     static constexpr char const *nameModifier = "Phase2";
 
@@ -483,6 +484,7 @@ namespace pixelTopology {
     static constexpr uint16_t lastRowInModule = numRowsInModule - 1;
     static constexpr uint16_t lastColInModule = numColsInModule - 1;
 
+    // 418 bins < 512, 9 bits are enough
     static constexpr uint16_t clusterBinning = numColsInModule + 2;
     static constexpr uint16_t clusterBits = 9;
 
@@ -557,7 +559,7 @@ namespace pixelTopology {
     static constexpr uint32_t maxPixInModule = 10000;
 
     static constexpr uint32_t maxNumOfActiveDoublets =
-        maxNumberOfDoublets / 4;  //TODO need to think a better way to avoid this duplication
+        maxNumberOfDoublets / 4;  // TODO need to think a better way to avoid this duplication
     static constexpr uint32_t maxCellsPerHit = 256;
 
     static constexpr uint32_t maxNumClustersPerModules = phase1HIonPixelTopology::maxNumClustersPerModules;

RecoLocalTracker/SiPixelClusterizer/plugins/SiPixelDigisClustersFromSoAAlpaka.cc

@@ -160,22 +160,26 @@ void SiPixelDigisClustersFromSoAAlpaka<TrackerTraits>::produce(edm::StreamID,
   std::cout << "Dumping all digis. nDigis = " << nDigis << std::endl;
 #endif
   for (uint32_t i = 0; i < nDigis; i++) {
-#ifdef GPU_DEBUG
-    PixelDigi dig2{digisView[i].pdigi()};
-    std::cout << i << ";" << digisView[i].rawIdArr() << ";" << digisView[i].clus() << ";" << digisView[i].pdigi() << ";"
-              << digisView[i].adc() << ";" << dig2.row() << ";" << dig2.column() << std::endl;
-#endif
-
     // check for uninitialized digis
     if (digisView[i].rawIdArr() == 0)
       continue;
     // check for noisy/dead pixels (electrons set to 0)
     if (digisView[i].adc() == 0)
       continue;
-    if (digisView[i].clus() >= -pixelClustering::invalidClusterId)
-      continue;  // not in cluster; TODO add an assert for the size
+    // not in cluster; TODO add an assert for the size
+    if (digisView[i].clus() == pixelClustering::invalidClusterId) {
+      continue;
+    }
+    // unexpected invalid value
+    if (digisView[i].clus() < pixelClustering::invalidClusterId) {
+      edm::LogError("SiPixelDigisClustersFromSoAAlpaka")
+          << "Skipping pixel digi with unexpected invalid cluster id " << digisView[i].clus();
+      continue;
+    }
+    // from clusters killed by charge cut
     if (digisView[i].clus() == pixelClustering::invalidModuleId)
-      continue;  // from clusters killed by charge cut
+      continue;
+
 #ifdef EDM_ML_DEBUG
     assert(digisView[i].rawIdArr() > 109999);
 #endif
@@ -200,6 +204,10 @@ void SiPixelDigisClustersFromSoAAlpaka<TrackerTraits>::produce(edm::StreamID,
       }
     }
     PixelDigi dig{digisView[i].pdigi()};
+#ifdef GPU_DEBUG
+    std::cout << i << ";" << digisView[i].rawIdArr() << ";" << digisView[i].clus() << ";" << digisView[i].pdigi() << ";"
+              << digisView[i].adc() << ";" << dig.row() << ";" << dig.column() << std::endl;
+#endif
 
     if (storeDigis_)
       (*detDigis).data.emplace_back(dig);

RecoLocalTracker/SiPixelClusterizer/plugins/SiPixelRawToClusterGPUKernel.cu

@@ -607,8 +607,8 @@ namespace pixelgpudetails {
           digis_d->moduleId(), clusters_d->moduleStart(), digis_d->clus(), wordCounter);
       cudaCheck(cudaGetLastError());
 
-      threadsPerBlock = ((TrackerTraits::maxPixInModule / 16 + 128 - 1) / 128) *
-                        128;  /// should be larger than maxPixInModule/16 aka (maxPixInModule/maxiter in the kernel)
+      // should be larger than maxPixInModule/16 aka (maxPixInModule/maxiter in the kernel)
+      threadsPerBlock = ((TrackerTraits::maxPixInModule / 16 + 128 - 1) / 128) * 128;
       blocks = TrackerTraits::numberOfModules;
 #ifdef GPU_DEBUG
       std::cout << "CUDA findClus kernel launch with " << blocks << " blocks of " << threadsPerBlock << " threads\n";
@@ -752,13 +752,13 @@ namespace pixelgpudetails {
     cudaCheck(cudaGetLastError());
 
     auto nModules_Clusters_d = cms::cuda::make_device_unique<uint32_t[]>(3, stream);
-    // MUST be ONE block
 
 #ifdef GPU_DEBUG
     cudaCheck(cudaStreamSynchronize(stream));
     std::cout << "CUDA fillHitsModuleStart kernel launch \n";
 #endif
 
+    // MUST be ONE block
     fillHitsModuleStart<TrackerTraits><<<1, 1024, 0, stream>>>(clusters_d->clusInModule(),
                                                                clusters_d->clusModuleStart(),
                                                                clusters_d->moduleStart(),

RecoLocalTracker/SiPixelClusterizer/plugins/alpaka/CalibPixel.h

@@ -23,9 +23,7 @@
 namespace calibPixel {
   using namespace cms::alpakatools;
 
-  constexpr uint16_t InvId = std::numeric_limits<uint16_t>::max() - 1;
-  // must be > MaxNumModules
-
+  template <bool debug = false>
   struct CalibDigis {
     template <typename TAcc>
     ALPAKA_FN_ACC void operator()(const TAcc& acc,
@@ -41,52 +39,57 @@ namespace calibPixel {
 
       // zero for next kernels...
       if (cms::alpakatools::once_per_grid(acc)) {
-        clus_view[0].clusModuleStart() = clus_view[0].moduleStart() = 0;
+        clus_view[0].clusModuleStart() = 0;
+        clus_view[0].moduleStart() = 0;
+      }
+      for (auto i : cms::alpakatools::elements_with_stride(acc, phase1PixelTopology::numberOfModules)) {
+        clus_view[i].clusInModule() = 0;
       }
 
-      cms::alpakatools::for_each_element_in_grid_strided(
-          acc, phase1PixelTopology::numberOfModules, [&](uint32_t i) { clus_view[i].clusInModule() = 0; });
-      cms::alpakatools::for_each_element_in_grid_strided(acc, numElements, [&](uint32_t i) {
+      for (auto i : cms::alpakatools::elements_with_stride(acc, numElements)) {
         auto dvgi = view[i];
-        if (dvgi.moduleId() != InvId) {
-          bool isDeadColumn = false, isNoisyColumn = false;
-          int row = dvgi.xx();
-          int col = dvgi.yy();
-          auto ret = SiPixelGainUtilities::getPedAndGain(gains, dvgi.moduleId(), col, row, isDeadColumn, isNoisyColumn);
-          float pedestal = ret.first;
-          float gain = ret.second;
-          if (isDeadColumn | isNoisyColumn) {
-            dvgi.moduleId() = InvId;
-            dvgi.adc() = 0;
+        if (dvgi.moduleId() == ::pixelClustering::invalidModuleId)
+          continue;
+
+        bool isDeadColumn = false, isNoisyColumn = false;
+        int row = dvgi.xx();
+        int col = dvgi.yy();
+        auto ret = SiPixelGainUtilities::getPedAndGain(gains, dvgi.moduleId(), col, row, isDeadColumn, isNoisyColumn);
+        float pedestal = ret.first;
+        float gain = ret.second;
+        if (isDeadColumn | isNoisyColumn) {
+          if constexpr (debug)
             printf("bad pixel at %d in %d\n", i, dvgi.moduleId());
-          } else {
-            float vcal = dvgi.adc() * gain - pedestal * gain;
+          dvgi.moduleId() = ::pixelClustering::invalidModuleId;
+          dvgi.adc() = 0;
+        } else {
+          float vcal = dvgi.adc() * gain - pedestal * gain;
 
-            float conversionFactor = dvgi.moduleId() < 96 ? VCaltoElectronGain_L1 : VCaltoElectronGain;
-            float offset = dvgi.moduleId() < 96 ? VCaltoElectronOffset_L1 : VCaltoElectronOffset;
+          float conversionFactor = dvgi.moduleId() < 96 ? VCaltoElectronGain_L1 : VCaltoElectronGain;
+          float offset = dvgi.moduleId() < 96 ? VCaltoElectronOffset_L1 : VCaltoElectronOffset;
 #ifdef GPU_DEBUG
-            auto old_adc = dvgi.adc();
+          auto old_adc = dvgi.adc();
 #endif
-            dvgi.adc() = std::max(100, int(vcal * conversionFactor + offset));
+          dvgi.adc() = std::max(100, int(vcal * conversionFactor + offset));
 #ifdef GPU_DEBUG
-            if (cms::alpakatools::once_per_grid(acc)) {
-              printf(
-                  "module %d pixel %d -> old_adc = %d; vcal = %.2f; conversionFactor = %.2f; offset = %.2f; new_adc = "
-                  "%d \n",
-                  dvgi.moduleId(),
-                  i,
-                  old_adc,
-                  vcal,
-                  conversionFactor,
-                  offset,
-                  dvgi.adc());
-            }
-#endif
+          if (cms::alpakatools::once_per_grid(acc)) {
+            printf(
+                "module %d pixel %d -> old_adc = %d; vcal = %.2f; conversionFactor = %.2f; offset = %.2f; new_adc = "
+                "%d \n",
+                dvgi.moduleId(),
+                i,
+                old_adc,
+                vcal,
+                conversionFactor,
+                offset,
+                dvgi.adc());
           }
+#endif
         }
-      });
+      }
     }
   };
+
   struct CalibDigisPhase2 {
     template <typename TAcc>
     ALPAKA_FN_ACC void operator()(const TAcc& acc,