Merge pull request #34882 from czangela/cpe-refactorings-2021-08

Follow up to PR review of CPEFast to better reproduce Generic and #34400
cms-sw · Aug 28, 2021 · 08f3fae · 08f3fae
2 parents 5c28861 + 6bc85e6
commit 08f3fae
Show file tree

Hide file tree

Showing 4 changed files with 57 additions and 34 deletions.
diff --git a/Geometry/TrackerGeometryBuilder/interface/phase1PixelTopology.h b/Geometry/TrackerGeometryBuilder/interface/phase1PixelTopology.h
@@ -19,6 +19,10 @@ namespace phase1PixelTopology {
   constexpr int16_t xOffset = -81;
   constexpr int16_t yOffset = -54 * 4;
 
+  constexpr uint16_t pixelThickness = 285;
+  constexpr uint16_t pixelPitchX = 100;
+  constexpr uint16_t pixelPitchY = 150;
+
   constexpr uint32_t numPixsInModule = uint32_t(numRowsInModule) * uint32_t(numColsInModule);
 
   constexpr uint32_t numberOfModules = 1856;

diff --git a/RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforGPU.h b/RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforGPU.h
@@ -12,6 +12,15 @@
 #include "HeterogeneousCore/CUDAUtilities/interface/cudaCompat.h"
 #include "CUDADataFormats/TrackingRecHit/interface/SiPixelHitStatus.h"
 
+namespace CPEFastParametrisation {
+  // From https://cmssdt.cern.ch/dxr/CMSSW/source/CondFormats/SiPixelTransient/src/SiPixelGenError.cc#485-486
+  // qbin: int (0-4) describing the charge of the cluster
+  // [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
+  constexpr int kGenErrorQBins = 5;
+  // arbitrary number of bins for sampling errors
+  constexpr int kNumErrorBins = 16;
+}  // namespace CPEFastParametrisation
+
 namespace pixelCPEforGPU {
 
   using Status = SiPixelHitStatus;
@@ -43,9 +52,10 @@ namespace pixelCPEforGPU {
 
     float apeXX, apeYY;  // ape^2
     uint8_t sx2, sy1, sy2;
-    uint8_t sigmax[16], sigmax1[16], sigmay[16];  // in micron
-    float xfact[5], yfact[5];
-    int minCh[5];
+    uint8_t sigmax[CPEFastParametrisation::kNumErrorBins], sigmax1[CPEFastParametrisation::kNumErrorBins],
+        sigmay[CPEFastParametrisation::kNumErrorBins];  // in micron
+    float xfact[CPEFastParametrisation::kGenErrorQBins], yfact[CPEFastParametrisation::kGenErrorQBins];
+    int minCh[CPEFastParametrisation::kGenErrorQBins];
 
     Frame frame;
   };
@@ -342,19 +352,26 @@ namespace pixelCPEforGPU {
 
     auto ch = cp.charge[ic];
     auto bin = 0;
-    for (; bin < 4; ++bin)
+    for (; bin < CPEFastParametrisation::kGenErrorQBins - 1; ++bin)
+      // find first bin which minimum charge exceeds cluster charge
       if (ch < detParams.minCh[bin + 1])
         break;
-    assert(bin < 5);
 
-    cp.status[ic].qBin = 4 - bin;
+    // in detParams qBins are reversed bin0 -> smallest charge, bin4-> largest charge
+    // whereas in CondFormats/SiPixelTransient/src/SiPixelGenError.cc it is the opposite
+    // so we reverse the bin here -> kGenErrorQBins - 1 - bin
+    cp.status[ic].qBin = CPEFastParametrisation::kGenErrorQBins - 1 - bin;
     cp.status[ic].isOneX = isOneX;
     cp.status[ic].isBigX = (isOneX & isBigX) | isEdgeX;
     cp.status[ic].isOneY = isOneY;
     cp.status[ic].isBigY = (isOneY & isBigY) | isEdgeY;
 
-    auto xoff = 81.f * comParams.thePitchX;
-    int jx = std::min(15, std::max(0, int(16.f * (cp.xpos[ic] + xoff) / (2 * xoff))));
+    auto xoff = -float(phase1PixelTopology::xOffset) * comParams.thePitchX;
+    int low_value = 0;
+    int high_value = CPEFastParametrisation::kNumErrorBins - 1;
+    int bin_value = float(CPEFastParametrisation::kNumErrorBins) * (cp.xpos[ic] + xoff) / (2 * xoff);
+    // return estimated bin value truncated to [0, 15]
+    int jx = std::clamp(bin_value, low_value, high_value);
 
     auto toCM = [](uint8_t x) { return float(x) * 1.e-4; };
 

diff --git a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFast.cc b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFast.cc
@@ -178,11 +178,14 @@ void PixelCPEFast::fillParamsForGpu() {
 
     auto toMicron = [&](float x) { return std::min(511, int(x * 1.e4f + 0.5f)); };
 
+    // average angle
+    auto gvx = p.theOrigin.x() + 40.f * commonParamsGPU_.thePitchX;
+    auto gvy = p.theOrigin.y();
+    auto gvz = 1.f / p.theOrigin.z();
+    //--- Note that the normalization is not required as only the ratio used
+
     {
-      // average angle
-      auto gvx = p.theOrigin.x() + 40.f * commonParamsGPU_.thePitchX;
-      auto gvy = p.theOrigin.y();
-      auto gvz = 1.f / p.theOrigin.z();
+      // calculate angles (fed into errorFromTemplates)
       cp.cotalpha = gvx * gvz;
       cp.cotbeta = gvy * gvz;
       errorFromTemplates(p, cp, 20000.);
@@ -204,9 +207,9 @@ void PixelCPEFast::fillParamsForGpu() {
     g.sy2 = std::max(55, toMicron(cp.sy2));  // sometimes sy2 is smaller than others (due to angle?)
 
     // sample xerr as function of position
-    auto const xoff = -81.f * commonParamsGPU_.thePitchX;
+    auto const xoff = float(phase1PixelTopology::xOffset) * commonParamsGPU_.thePitchX;
 
-    for (int ix = 0; ix < 16; ++ix) {
+    for (int ix = 0; ix < CPEFastParametrisation::kNumErrorBins; ++ix) {
       auto x = xoff * (1.f - (0.5f + float(ix)) / 8.f);
       auto gvx = p.theOrigin.x() - x;
       auto gvy = p.theOrigin.y();
@@ -221,8 +224,8 @@ void PixelCPEFast::fillParamsForGpu() {
     }
 #ifdef EDM_ML_DEBUG
     // sample yerr as function of position
-    auto const yoff = -54.f * 4.f * commonParamsGPU_.thePitchY;
-    for (int ix = 0; ix < 16; ++ix) {
+    auto const yoff = float(phase1PixelTopology::yOffset) * commonParamsGPU_.thePitchY;
+    for (int ix = 0; ix < CPEFastParametrisation::kNumErrorBins; ++ix) {
       auto y = yoff * (1.f - (0.5f + float(ix)) / 8.f);
       auto gvx = p.theOrigin.x() + 40.f * commonParamsGPU_.thePitchY;
       auto gvy = p.theOrigin.y() - y;
@@ -235,19 +238,13 @@ void PixelCPEFast::fillParamsForGpu() {
     }
 #endif  // EDM_ML_DEBUG
 
-    // average angle
-    auto gvx = p.theOrigin.x() + 40.f * commonParamsGPU_.thePitchX;
-    auto gvy = p.theOrigin.y();
-    auto gvz = 1.f / p.theOrigin.z();
-    //--- Note that the normalization is not required as only the ratio used
-
-    // calculate angles (fed into errorFromTemplates)
+    // calculate angles (repeated)
     cp.cotalpha = gvx * gvz;
     cp.cotbeta = gvy * gvz;
     auto aveCB = cp.cotbeta;
 
     // sample x by charge
-    int qbin = 5;  // low charge
+    int qbin = CPEFastParametrisation::kGenErrorQBins;  // low charge
     int k = 0;
     for (int qclus = 1000; qclus < 200000; qclus += 1000) {
       errorFromTemplates(p, cp, qclus);
@@ -263,25 +260,26 @@ void PixelCPEFast::fillParamsForGpu() {
                                << m * cp.sigmay << ' ' << m * cp.sy1 << ' ' << m * cp.sy2 << std::endl;
 #endif  // EDM_ML_DEBUG
     }
-    assert(k <= 5);
+    assert(k <= CPEFastParametrisation::kGenErrorQBins);
     // fill the rest  (sometimes bin 4 is missing)
-    for (int kk = k; kk < 5; ++kk) {
+    for (int kk = k; kk < CPEFastParametrisation::kGenErrorQBins; ++kk) {
       g.xfact[kk] = g.xfact[k - 1];
       g.yfact[kk] = g.yfact[k - 1];
       g.minCh[kk] = g.minCh[k - 1];
     }
     auto detx = 1.f / g.xfact[0];
     auto dety = 1.f / g.yfact[0];
-    for (int kk = 0; kk < 5; ++kk) {
+    for (int kk = 0; kk < CPEFastParametrisation::kGenErrorQBins; ++kk) {
       g.xfact[kk] *= detx;
       g.yfact[kk] *= dety;
     }
     // sample y in "angle"  (estimated from cluster size)
     float ys = 8.f - 4.f;  // apperent bias of half pixel (see plot)
+    // plot: https://indico.cern.ch/event/934821/contributions/3974619/attachments/2091853/3515041/DigilessReco.pdf page 25
     // sample yerr as function of "size"
-    for (int iy = 0; iy < 16; ++iy) {
+    for (int iy = 0; iy < CPEFastParametrisation::kNumErrorBins; ++iy) {
       ys += 1.f;  // first bin 0 is for size 9  (and size is in fixed point 2^3)
-      if (15 == iy)
+      if (CPEFastParametrisation::kNumErrorBins - 1 == iy)
         ys += 8.f;  // last bin for "overflow"
       // cp.cotalpha = ys*(commonParamsGPU_.thePitchX/(8.f*thickness));  //  use this to print sampling in "x"  (and comment the line below)
       cp.cotbeta = std::copysign(ys * (commonParamsGPU_.thePitchY / (8.f * thickness)), aveCB);

diff --git a/RecoPixelVertexing/PixelTriplets/plugins/BrokenLineFitOnGPU.h b/RecoPixelVertexing/PixelTriplets/plugins/BrokenLineFitOnGPU.h
@@ -90,16 +90,20 @@ __global__ void kernel_BLFastFit(Tuples const *__restrict__ foundNtuplets,
 #ifdef YERR_FROM_DC
       auto const &dp = hhp->cpeParams().detParams(hhp->detectorIndex(hit));
       auto status = hhp->status(hit);
-      int qbin = 4 - status.qBin;
+      int qbin = CPEFastParametrisation::kGenErrorQBins - 1 - status.qBin;
       assert(qbin >= 0 && qbin < 5);
       bool nok = (status.isBigY | status.isOneY);
       // compute cotanbeta and use it to recompute error
       dp.frame.rotation().multiply(dx, dy, dz, ux, uy, uz);
       auto cb = std::abs(uy / uz);
-      int bin = int(cb * (285.f / 150.f) * 8.f) - 4;
-      bin = std::max(0, std::min(15, bin));
-      float yerr = dp.sigmay[bin] * 1.e-4f;
-      yerr *= dp.yfact[qbin];  // inflate
+      int bin =
+          int(cb * (float(phase1PixelTopology::pixelThickess) / float(phase1PixelTopology::pixelPitchY)) * 8.f) - 4;
+      int low_value = 0;
+      int high_value = CPEFastParametrisation::kNumErrorBins - 1;
+      // return estimated bin value truncated to [0, 15]
+      bin = std::clamp(bin, low_value, high_value);
+      float yerr = dp.sigmay[bin] * 1.e-4f;  // toCM
+      yerr *= dp.yfact[qbin];                // inflate
       yerr *= yerr;
       yerr += dp.apeYY;
       yerr = nok ? hhp->yerrLocal(hit) : yerr;