Phase II Patatrack Pixel Local Reco

Author: AdrianoDee

CUDADataFormats/SiPixelCluster/interface/gpuClusteringConstants.h

@@ -11,11 +11,13 @@ namespace gpuClustering {
 #else
   // optimized for real data PU 50
   // tested on MC events with 55-75 pileup events
-  constexpr uint32_t maxHitsInIter() { return 160; }
+  constexpr uint32_t maxHitsInIter() { return 160; } //TODO better tuning for PU 140-200
 #endif
   constexpr uint32_t maxHitsInModule() { return 1024; }
 
-  constexpr uint16_t maxNumModules = 2000;
+  constexpr uint32_t maxNumDigis = 3*256*1024; // @PU=200 µ=530k σ=50k this is >4σ away
+  constexpr uint16_t maxNumModules = 4000;
+
   constexpr int32_t maxNumClustersPerModules = maxHitsInModule();
   constexpr uint16_t invalidModuleId = std::numeric_limits<uint16_t>::max() - 1;
   constexpr int invalidClusterId = -9999;

CUDADataFormats/TrackingRecHit/interface/TrackingRecHit2DHeterogeneous.h

@@ -3,6 +3,7 @@
 
 #include "CUDADataFormats/TrackingRecHit/interface/TrackingRecHit2DSOAView.h"
 #include "CUDADataFormats/Common/interface/HeterogeneousSoA.h"
+#include "RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforGPU.h"
 
 template <typename Traits>
 class TrackingRecHit2DHeterogeneous {
@@ -16,6 +17,7 @@ class TrackingRecHit2DHeterogeneous {
 
   explicit TrackingRecHit2DHeterogeneous(
       uint32_t nHits,
+      bool isUpgrade,
       int32_t offsetBPIX2,
       pixelCPEforGPU::ParamsOnGPU const* cpeParams,
       uint32_t const* hitsModuleStart,
@@ -33,6 +35,7 @@ class TrackingRecHit2DHeterogeneous {
   TrackingRecHit2DSOAView const* view() const { return m_view.get(); }
 
   auto nHits() const { return m_nHits; }
+  auto nMaxModules() const { return m_nMaxModules; }
   auto offsetBPIX2() const { return m_offsetBPIX2; }
 
   auto hitsModuleStart() const { return m_hitsModuleStart; }
@@ -62,10 +65,12 @@ class TrackingRecHit2DHeterogeneous {
   unique_ptr<TrackingRecHit2DSOAView> m_view;  //!
 
   uint32_t m_nHits;
+  bool m_isUpgrade;
   int32_t m_offsetBPIX2;
 
   uint32_t const* m_hitsModuleStart;  // needed for legacy, this is on GPU!
 
+  uint32_t m_nMaxModules;
   // needed as kernel params...
   PhiBinner* m_phiBinner;
   PhiBinner::index_type* m_phiBinnerStorage;
@@ -83,15 +88,18 @@ using TrackingRecHit2DHost = TrackingRecHit2DHeterogeneous<cms::cudacompat::Host
 template <typename Traits>
 TrackingRecHit2DHeterogeneous<Traits>::TrackingRecHit2DHeterogeneous(
     uint32_t nHits,
+    bool isUpgrade,
     int32_t offsetBPIX2,
     pixelCPEforGPU::ParamsOnGPU const* cpeParams,
     uint32_t const* hitsModuleStart,
     cudaStream_t stream,
     TrackingRecHit2DHeterogeneous<cms::cudacompat::GPUTraits> const* input)
-    : m_nHits(nHits), m_offsetBPIX2(offsetBPIX2), m_hitsModuleStart(hitsModuleStart) {
+    : m_nHits(nHits), m_isUpgrade(isUpgrade), m_offsetBPIX2(offsetBPIX2), m_hitsModuleStart(hitsModuleStart) {
   auto view = Traits::template make_host_unique<TrackingRecHit2DSOAView>(stream);
 
+  m_nMaxModules = m_isUpgrade ? phase2PixelTopology::numberOfModules : phase1PixelTopology::numberOfModules;
   view->m_nHits = nHits;
+  view->m_nMaxModules = m_nMaxModules;
   m_view = Traits::template make_unique<TrackingRecHit2DSOAView>(stream);  // leave it on host and pass it by value?
   m_AverageGeometryStore = Traits::template make_unique<TrackingRecHit2DSOAView::AverageGeometry>(stream);
   view->m_averageGeometry = m_AverageGeometryStore.get();
@@ -120,8 +128,9 @@ TrackingRecHit2DHeterogeneous<Traits>::TrackingRecHit2DHeterogeneous(
     copyFromGPU(input, stream);
   } else {
     assert(input == nullptr);
+    auto nL = m_isUpgrade ? phase2PixelTopology::numberOfLayers : phase1PixelTopology::numberOfLayers;
     m_store16 = Traits::template make_unique<uint16_t[]>(nHits * n16, stream);
-    m_store32 = Traits::template make_unique<float[]>(nHits * n32 + phase1PixelTopology::numberOfLayers + 1, stream);
+    m_store32 = Traits::template make_unique<float[]>(nHits * n32 + nL + 1, stream);
     m_PhiBinnerStore = Traits::template make_unique<TrackingRecHit2DSOAView::PhiBinner>(stream);
   }
 

CUDADataFormats/TrackingRecHit/interface/TrackingRecHit2DSOAView.h

@@ -6,7 +6,7 @@
 #include "CUDADataFormats/SiPixelCluster/interface/gpuClusteringConstants.h"
 #include "HeterogeneousCore/CUDAUtilities/interface/HistoContainer.h"
 #include "HeterogeneousCore/CUDAUtilities/interface/cudaCompat.h"
-#include "Geometry/TrackerGeometryBuilder/interface/phase1PixelTopology.h"
+#include "Geometry/TrackerGeometryBuilder/interface/pixelTopology.h"
 #include "CUDADataFormats/TrackingRecHit/interface/SiPixelHitStatus.h"
 
 namespace pixelCPEforGPU {
@@ -20,15 +20,17 @@ class TrackingRecHit2DSOAView {
 
   using hindex_type = uint32_t;  // if above is <=2^32
 
-  using PhiBinner = cms::cuda::HistoContainer<int16_t, 128, -1, 8 * sizeof(int16_t), hindex_type, 10>;
+  using PhiBinner =
+      cms::cuda::HistoContainer<int16_t, 256, -1, 8 * sizeof(int16_t), hindex_type, 28>;  //28 for phase2 geometry
 
-  using AverageGeometry = phase1PixelTopology::AverageGeometry;
+  using AverageGeometry = pixelTopology::AverageGeometry;
 
   template <typename>
   friend class TrackingRecHit2DHeterogeneous;
   friend class TrackingRecHit2DReduced;
 
   __device__ __forceinline__ uint32_t nHits() const { return m_nHits; }
+  __device__ __forceinline__ uint32_t nMaxModules() const { return m_nMaxModules; }
 
   __device__ __forceinline__ float& xLocal(int i) { return m_xl[i]; }
   __device__ __forceinline__ float xLocal(int i) const { return __ldg(m_xl + i); }
@@ -114,6 +116,7 @@ class TrackingRecHit2DSOAView {
   PhiBinner::index_type* m_phiBinnerStorage;
 
   uint32_t m_nHits;
+  uint32_t m_nMaxModules;
 };
 
 #endif  // CUDADataFormats_TrackingRecHit_interface_TrackingRecHit2DSOAView_h

CUDADataFormats/TrackingRecHit/src/TrackingRecHit2DHeterogeneous.cc

@@ -13,9 +13,9 @@ cms::cuda::host::unique_ptr<float[]> TrackingRecHit2DGPU::localCoordToHostAsync(
 
 template <>
 cms::cuda::host::unique_ptr<uint32_t[]> TrackingRecHit2DGPU::hitsModuleStartToHostAsync(cudaStream_t stream) const {
-  auto ret = cms::cuda::make_host_unique<uint32_t[]>(gpuClustering::maxNumModules + 1, stream);
+  auto ret = cms::cuda::make_host_unique<uint32_t[]>(nMaxModules() + 1, stream);
   cudaCheck(cudaMemcpyAsync(
-      ret.get(), m_hitsModuleStart, sizeof(uint32_t) * (gpuClustering::maxNumModules + 1), cudaMemcpyDefault, stream));
+      ret.get(), m_hitsModuleStart, sizeof(uint32_t) * (nMaxModules() + 1), cudaMemcpyDefault, stream));
   return ret;
 }
 

CUDADataFormats/TrackingRecHit/test/TrackingRecHit2DCUDA_t.cpp

@@ -16,16 +16,18 @@ int main() {
   cudaCheck(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking));
 
   auto nHits = 200;
+  auto nModules = 2000;
   // inner scope to deallocate memory before destroying the stream
   {
-    TrackingRecHit2DGPU tkhit(nHits, 0, nullptr, nullptr, stream);
+    TrackingRecHit2DGPU tkhit(nHits, nModules, 0, nullptr, nullptr, stream);
 
     testTrackingRecHit2D::runKernels(tkhit.view());
 
-    TrackingRecHit2DHost tkhitH(nHits, 0, nullptr, nullptr, stream, &tkhit);
+    TrackingRecHit2DHost tkhitH(nHits, nModules, 0, nullptr, nullptr, stream, &tkhit);
     cudaStreamSynchronize(stream);
     assert(tkhitH.view());
     assert(tkhitH.view()->nHits() == unsigned(nHits));
+    assert(tkhitH.view()->nMaxModules() == unsigned(nModules));
   }
 
   cudaCheck(cudaStreamDestroy(stream));

CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLWorker.cc

@@ -45,7 +45,7 @@
 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyMap.h"
 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
-#include "Geometry/TrackerGeometryBuilder/interface/phase1PixelTopology.h"
+#include "Geometry/TrackerGeometryBuilder/interface/pixelTopology.h"
 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"

…yBuilder/interface/phase1PixelTopology.h …eometryBuilder/interface/pixelTopology.hGeometry/TrackerGeometryBuilder/interface/phase1PixelTopology.h renamed to Geometry/TrackerGeometryBuilder/interface/pixelTopology.h Geometry/TrackerGeometryBuilder/interface/phase1PixelTopology.h renamed to Geometry/TrackerGeometryBuilder/interface/pixelTopology.h

@@ -1,11 +1,52 @@
-#ifndef Geometry_TrackerGeometryBuilder_phase1PixelTopology_h
-#define Geometry_TrackerGeometryBuilder_phase1PixelTopology_h
+#ifndef Geometry_TrackerGeometryBuilder_pixelTopology_h
+#define Geometry_TrackerGeometryBuilder_pixelTopology_h
 
 #include <cstdint>
 #include <array>
 
+namespace pixelTopology {
+  template <class Function, std::size_t... Indices>
+  constexpr auto map_to_array_helper(Function f, std::index_sequence<Indices...>)
+      -> std::array<typename std::result_of<Function(std::size_t)>::type, sizeof...(Indices)> {
+    return {{f(Indices)...}};
+  }
+
+  template <int N, class Function>
+  constexpr auto map_to_array(Function f) -> std::array<typename std::result_of<Function(std::size_t)>::type, N> {
+    return map_to_array_helper(f, std::make_index_sequence<N>{});
+  }
+
+  static constexpr auto maxNumberOfLadders = 160;
+  constexpr uint32_t maxLayers = 28;
+
+  struct AverageGeometry {
+    //
+    float ladderZ[maxNumberOfLadders];
+    float ladderX[maxNumberOfLadders];
+    float ladderY[maxNumberOfLadders];
+    float ladderR[maxNumberOfLadders];
+    float ladderMinZ[maxNumberOfLadders];
+    float ladderMaxZ[maxNumberOfLadders];
+    float endCapZ[2];  // just for pos and neg Layer1
+  };
+
+  constexpr inline uint16_t localY(uint16_t py, uint16_t n) {
+    auto roc = py / n;
+    auto shift = 2 * roc;
+    auto yInRoc = py - n * roc;
+    if (yInRoc > 0)
+      shift += 1;
+    return py + shift;
+  }
+
+}  // namespace pixelTopology
+
 namespace phase1PixelTopology {
 
+  constexpr uint16_t numberOfModulesInBarrel = 1184;
+  constexpr uint16_t numberOfModulesInLadder = 8;
+  constexpr uint16_t numberOfLaddersInBarrel = numberOfModulesInBarrel / numberOfModulesInLadder;
+
   constexpr uint16_t numRowsInRoc = 80;
   constexpr uint16_t numColsInRoc = 52;
   constexpr uint16_t lastRowInRoc = numRowsInRoc - 1;
@@ -54,21 +95,7 @@ namespace phase1PixelTopology {
       "E-3"  // negative endcap
   };
 
-  constexpr uint32_t numberOfModulesInBarrel = 1184;
-  constexpr uint32_t numberOfLaddersInBarrel = numberOfModulesInBarrel / 8;
-
-  template <class Function, std::size_t... Indices>
-  constexpr auto map_to_array_helper(Function f, std::index_sequence<Indices...>)
-      -> std::array<typename std::result_of<Function(std::size_t)>::type, sizeof...(Indices)> {
-    return {{f(Indices)...}};
-  }
-
-  template <int N, class Function>
-  constexpr auto map_to_array(Function f) -> std::array<typename std::result_of<Function(std::size_t)>::type, N> {
-    return map_to_array_helper(f, std::make_index_sequence<N>{});
-  }
-
-  constexpr uint32_t findMaxModuleStride() {
+  constexpr uint16_t findMaxModuleStride() {
     bool go = true;
     int n = 2;
     while (go) {
@@ -85,7 +112,7 @@ namespace phase1PixelTopology {
     return n / 2;
   }
 
-  constexpr uint32_t maxModuleStride = findMaxModuleStride();
+  constexpr uint16_t maxModuleStride = findMaxModuleStride();
 
   constexpr uint8_t findLayer(uint32_t detId, uint8_t sl = 0) {
     for (uint8_t i = sl; i < std::size(layerStart); ++i)
@@ -107,7 +134,7 @@ namespace phase1PixelTopology {
   __device__
 #endif
       constexpr std::array<uint8_t, layerIndexSize>
-          layer = map_to_array<layerIndexSize>(findLayerFromCompact);
+          layer = pixelTopology::map_to_array<layerIndexSize>(findLayerFromCompact); 
 
   constexpr uint8_t getLayer(uint32_t detId) {
     return phase1PixelTopology::layer[detId / phase1PixelTopology::maxModuleStride];
@@ -117,7 +144,7 @@ namespace phase1PixelTopology {
     bool res = true;
     for (auto i = 0U; i < numberOfModules; ++i) {
       auto j = i / maxModuleStride;
-      res &= (layer[j] < 10);
+      res &= (layer[j] < numberOfLayers);
       res &= (i >= layerStart[layer[j]]);
       res &= (i < layerStart[layer[j] + 1]);
     }
@@ -172,18 +199,67 @@ namespace phase1PixelTopology {
     return py + shift;
   }
 
-  //FIXME move it elsewhere?
-  struct AverageGeometry {
-    static constexpr auto numberOfLaddersInBarrel = phase1PixelTopology::numberOfLaddersInBarrel;
-    float ladderZ[numberOfLaddersInBarrel];
-    float ladderX[numberOfLaddersInBarrel];
-    float ladderY[numberOfLaddersInBarrel];
-    float ladderR[numberOfLaddersInBarrel];
-    float ladderMinZ[numberOfLaddersInBarrel];
-    float ladderMaxZ[numberOfLaddersInBarrel];
-    float endCapZ[2];  // just for pos and neg Layer1
-  };
-
 }  // namespace phase1PixelTopology
 
-#endif  // Geometry_TrackerGeometryBuilder_phase1PixelTopology_h
+namespace phase2PixelTopology {
+
+  constexpr uint32_t numberOfModulesInBarrel = 756;
+  constexpr uint32_t numberOfModulesInLadder = 9;
+  constexpr uint32_t numberOfLaddersInBarrel = numberOfModulesInBarrel / numberOfModulesInLadder;
+
+  constexpr uint32_t numberOfModules = 3892;
+  constexpr uint8_t numberOfLayers = 28;
+
+  constexpr uint32_t layerStart[numberOfLayers + 1] = {
+      0,    108,  324,  504,  //Barrel
+      756,  864,  972,  1080, 1188, 1296, 1404, 1512, 1620, 1796, 1972, 2148, //Fp
+      2324, 2432, 2540, 2648, 2756, 2864, 2972, 3080, 3188, 3364, 3540, 3716, //Np
+      numberOfModules};
+
+  constexpr uint16_t findMaxModuleStride() {
+    bool go = true;
+    int n = 2;
+    while (go) {
+      for (uint8_t i = 1; i < numberOfLayers + 1; ++i) {
+        if (layerStart[i] % n != 0) {
+          go = false;
+          break;
+        }
+      }
+      if (!go)
+        break;
+      n *= 2;
+    }
+    return n / 2;
+  }
+
+  constexpr uint16_t maxModuleStride = findMaxModuleStride();
+
+  constexpr uint8_t findLayerFromCompact(uint32_t detId) {
+    detId *= maxModuleStride;
+    for (uint8_t i = 0; i < numberOfLayers + 1; ++i)
+      if (detId < layerStart[i + 1])
+        return i;
+    return numberOfLayers + 1;
+  }
+
+  constexpr uint16_t layerIndexSize = numberOfModules / maxModuleStride;
+  constexpr std::array<uint8_t, layerIndexSize> layer =
+      pixelTopology::map_to_array<layerIndexSize>(findLayerFromCompact);
+
+  constexpr bool validateLayerIndex() {
+    bool res = true;
+    for (auto i = 0U; i < numberOfModules; ++i) {
+      auto j = i / maxModuleStride;
+      res &= (layer[j] < numberOfLayers);
+      res &= (i >= layerStart[layer[j]]);
+      res &= (i < layerStart[layer[j] + 1]);
+    }
+    return res;
+  }
+
+  static_assert(validateLayerIndex(), "phase2 layer from detIndex algo is buggy");
+
+}  // namespace phase2PixelTopology
+
+#endif  // Geometry_TrackerGeometryBuilder_pixelTopology_h

Geometry/TrackerGeometryBuilder/test/phase1PixelTopology_t.cu

@@ -2,7 +2,7 @@
 #include <iostream>
 #include <tuple>
 
-#include "Geometry/TrackerGeometryBuilder/interface/phase1PixelTopology.h"
+#include "Geometry/TrackerGeometryBuilder/interface/pixelTopology.h"
 
 #include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h"
 #include "HeterogeneousCore/CUDAUtilities/interface/requireDevices.h"

RecoLocalTracker/SiPixelClusterizer/plugins/SiPixelClusterThresholds.h

@@ -10,5 +10,6 @@ struct SiPixelClusterThresholds {
 };
 
 constexpr SiPixelClusterThresholds kSiPixelClusterThresholdsDefaultPhase1{.layer1 = 2000, .otherLayers = 4000};
+constexpr SiPixelClusterThresholds kSiPixelClusterThresholdsDefaultPhase2{.layer1 = 4000, .otherLayers = 4000};
 
 #endif  // RecoLocalTracker_SiPixelClusterizer_plugins_SiPixelClusterThresholds_h