Pixel cluster repair by morphing

Configuration/ProcessModifiers/python/siPixelDigiMorphing_cff.py

@@ -0,0 +1,4 @@
+import FWCore.ParameterSet.Config as cms
+
+# This modifier is to run the SiPixelDigiMorphing
+siPixelDigiMorphing =  cms.Modifier()

RecoLocalTracker/Configuration/python/RecoLocalTracker_cff.py

@@ -8,13 +8,17 @@
 from RecoLocalTracker.SiStripZeroSuppression.SiStripZeroSuppression_cfi import *
 from RecoLocalTracker.SiStripClusterizer.SiStripClusterizer_cfi import *
 from RecoLocalTracker.SiPixelClusterizer.siPixelClustersPreSplitting_cff import *
+from RecoLocalTracker.SiPixelDigiReProducers.siPixelDigisMorphed_cfi import *
 from RecoLocalTracker.SiPixelRecHits.SiPixelRecHits_cfi import *
 from RecoLocalTracker.SubCollectionProducers.clustersummaryproducer_cfi import *
 
 pixeltrackerlocalrecoTask = cms.Task(
     siPixelClustersPreSplittingTask,
     siPixelRecHitsPreSplittingTask)
 
+from Configuration.ProcessModifiers.siPixelDigiMorphing_cff import *
+siPixelDigiMorphing.toModify(pixeltrackerlocalrecoTask, func=lambda t: t.add(siPixelDigisMorphed))
+
 striptrackerlocalrecoTask = cms.Task(
     siStripZeroSuppression,
     siStripClusters,

RecoLocalTracker/SiPixelClusterizer/plugins/PixelThresholdClusterizer.cc

@@ -62,6 +62,7 @@ PixelThresholdClusterizer::PixelThresholdClusterizer(edm::ParameterSet const& co
       doMissCalibrate(conf.getParameter<bool>("MissCalibrate")),
       doSplitClusters(conf.getParameter<bool>("SplitClusters")) {
   theBuffer.setSize(theNumOfRows, theNumOfCols);
+  theFakePixels.clear();
 }
 /////////////////////////////////////////////////////////////////////////////
 PixelThresholdClusterizer::~PixelThresholdClusterizer() {}
@@ -110,6 +111,8 @@ bool PixelThresholdClusterizer::setup(const PixelGeomDetUnit* pixDet) {
     theBuffer.setSize(nrows, ncols);  // Modify
   }
 
+  theFakePixels.resize(nrows * ncols, false);
+
   return true;
 }
 //----------------------------------------------------------------------------
@@ -179,6 +182,8 @@ void PixelThresholdClusterizer::clusterizeDetUnitT(const T& input,
 
   //  Need to clean unused pixels from the buffer array.
   clear_buffer(begin, end);
+
+  theFakePixels.clear();
 }
 
 //----------------------------------------------------------------------------
@@ -218,7 +223,7 @@ void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator
   }
 #endif
   int electron[end - begin];  // pixel charge in electrons
-  memset(electron, 0, sizeof(electron));
+  memset(electron, 0, (end - begin) * sizeof(int));
 
   if (doPhase2Calibration) {
     int i = 0;
@@ -257,7 +262,9 @@ void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator
   for (DigiIterator di = begin; di != end; ++di) {
     int row = di->row();
     int col = di->column();
-    int adc = electron[i++];  // this is in electrons
+    // VV: do not calibrate a fake pixel, it already has a unit of 10e-:
+    int adc = (di->flag() != 0) ? di->adc() * 10 : electron[i];  // this is in electrons
+    i++;
 
 #ifdef PIXELREGRESSION
     int adcOld = calibrate(di->adc(), col, row);
@@ -277,6 +284,9 @@ void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator
 
     if (adc >= thePixelThreshold) {
       theBuffer.set_adc(row, col, adc);
+      // VV: add pixel to the fake list. Only when running on digi collection
+      if (di->flag() != 0)
+        theFakePixels[row * theNumOfCols + col] = true;
       if (adc >= theSeedThreshold)
         theSeeds.push_back(SiPixelCluster::PixelPos(row, col));
     }
@@ -414,8 +424,9 @@ SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::Pix
   theBuffer.set_adc(pix, 1);
   //  }
 
-  AccretionCluster acluster;
+  AccretionCluster acluster, cldata;
   acluster.add(pix, seed_adc);
+  cldata.add(pix, seed_adc);
 
   //Here we search all pixels adjacent to all pixels in the cluster.
   bool dead_flag = false;
@@ -433,6 +444,10 @@ SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::Pix
           SiPixelCluster::PixelPos newpix(r, c);
           if (!acluster.add(newpix, theBuffer(r, c)))
             goto endClus;
+          // VV: no fake pixels in cluster, leads to non-contiguous clusters
+          if (!theFakePixels[r * theNumOfCols + c]) {
+            cldata.add(newpix, theBuffer(r, c));
+          }
           theBuffer.set_adc(newpix, 1);
         }
 
@@ -464,7 +479,7 @@ SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::Pix
 
   }  // while accretion
 endClus:
-  SiPixelCluster cluster(acluster.isize, acluster.adc, acluster.x, acluster.y, acluster.xmin, acluster.ymin);
+  SiPixelCluster cluster(cldata.isize, cldata.adc, cldata.x, cldata.y, cldata.xmin, cldata.ymin);
   //Here we split the cluster, if the flag to do so is set and we have found a dead or noisy pixel.
 
   if (dead_flag && doSplitClusters) {

RecoLocalTracker/SiPixelClusterizer/plugins/PixelThresholdClusterizer.h

@@ -90,6 +90,8 @@ class dso_hidden PixelThresholdClusterizer final : public PixelClusterizerBase {
   std::vector<SiPixelCluster::PixelPos> theSeeds;  // cached seed pixels
   std::vector<SiPixelCluster> theClusters;         // resulting clusters
 
+  std::vector<bool> theFakePixels;  // fake pixels introduced to guide clustering
+
   //! Clustering-related quantities:
   float thePixelThresholdInNoiseUnits;    // Pixel threshold in units of noise
   float theSeedThresholdInNoiseUnits;     // Pixel cluster seed in units of noise

RecoLocalTracker/SiPixelClusterizer/python/SiPixelClusterizerPreSplitting_cfi.py

@@ -23,3 +23,11 @@
         )
     )
 )
+
+from Configuration.ProcessModifiers.siPixelDigiMorphing_cff import siPixelDigiMorphing
+siPixelDigiMorphing.toModify(
+    siPixelClustersPreSplitting,
+    cpu = dict(
+         src = 'siPixelDigisMorphed'
+    )
+)

RecoLocalTracker/SiPixelDigiReProducers/doc/SiPixelDigiReProducers.doc

@@ -0,0 +1,6 @@
+
+/*!
+
+\page RecoLocalTracker_SiPixelDigiReProducers Package RecoLocalTracker/SiPixelDigiReProducers
+*/
+

RecoLocalTracker/SiPixelDigiReProducers/plugins/BuildFile.xml

@@ -0,0 +1,7 @@
+<use name="DataFormats/Common"/>
+<use name="FWCore/ParameterSet"/>
+<use name="DataFormats/SiPixelDigi"/>
+<use name="CalibTracker/SiPixelESProducers"/>
+<library file="*.cc" name="RecoLocalTrackerSiPixelDigiReProducersPlugins">
+  <flags EDM_PLUGIN="1"/>
+</library>

RecoLocalTracker/SiPixelDigiReProducers/plugins/SiPixelDigiMorphing.cc

@@ -0,0 +1,207 @@
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "DataFormats/Common/interface/DetSetVector.h"
+#include "DataFormats/Common/interface/DetSetVectorNew.h"
+#include "DataFormats/DetId/interface/DetId.h"
+#include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
+
+#include <bitset>
+#include <iterator>
+
+class SiPixelDigiMorphing : public edm::stream::EDProducer<> {
+public:
+  explicit SiPixelDigiMorphing(const edm::ParameterSet& conf);
+  ~SiPixelDigiMorphing() override = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+  void produce(edm::Event& e, const edm::EventSetup& c) override;
+
+private:
+  edm::EDGetTokenT<edm::DetSetVector<PixelDigi>> tPixelDigi_;
+  edm::EDPutTokenT<edm::DetSetVector<PixelDigi>> tPutPixelDigi_;
+
+  const int32_t nrows_;
+  const int32_t ncols_;
+  const int32_t nrocs_;  // in Phase 1, this is ROCs, but could be any subset of a pixel row
+  const int32_t iters_;
+  const uint32_t fakeAdc_;
+
+  int32_t ncols_r_;  // number of columns per ROC
+  int32_t ksize_;    // kernel size
+
+  std::vector<uint64_t> kernel1_;
+  std::vector<uint64_t> kernel2_;
+  uint64_t mask_;
+
+  enum MorphOption { kDilate, kErode };
+
+  void morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const;
+};
+
+SiPixelDigiMorphing::SiPixelDigiMorphing(edm::ParameterSet const& conf)
+    : tPixelDigi_(consumes(conf.getParameter<edm::InputTag>("src"))),
+      tPutPixelDigi_(produces<edm::DetSetVector<PixelDigi>>()),
+      nrows_(conf.getParameter<int32_t>("nrows")),
+      ncols_(conf.getParameter<int32_t>("ncols")),
+      nrocs_(conf.getParameter<int32_t>("nrocs")),
+      iters_(conf.getParameter<int32_t>("iters")),
+      fakeAdc_(conf.getParameter<uint32_t>("fakeAdc")) {
+  if (ncols_ % nrocs_ == 0) {
+    ncols_r_ = ncols_ / nrocs_;
+  } else {
+    throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
+                                          << " number of columns not divisible with"
+                                          << " number of ROCs\n";
+  }
+
+  if (ncols_r_ + 2 * iters_ <= int(sizeof(uint64_t) * 8)) {
+    ksize_ = 2 * iters_ + 1;
+  } else {
+    throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
+                                          << " too many columns per ROC"
+                                          << " or too many iterations set\n"
+                                          << " Ncol/Nrocs+2*iters should not be"
+                                          << " more than " << sizeof(uint64_t) * 8 << "\n";
+  }
+
+  std::vector<int32_t> k1(conf.getParameter<std::vector<int32_t>>("kernel1"));
+  std::vector<int32_t> k2(conf.getParameter<std::vector<int32_t>>("kernel2"));
+
+  kernel1_.resize(ksize_, 0);
+  kernel2_.resize(ksize_, 0);
+  mask_ = 0;
+  int w = (ncols_r_ + 2 * iters_) / ksize_ + ((ncols_r_ + 2 * iters_) % ksize_ != 0);
+  for (int j = 0; j < w; j++) {
+    for (int ii = 0; ii < ksize_; ii++) {
+      kernel1_[ii] <<= ksize_;
+      kernel1_[ii] |= k1[ii];
+      kernel2_[ii] <<= ksize_;
+      kernel2_[ii] |= k2[ii];
+    }
+    mask_ <<= ksize_;
+    mask_ |= 1;
+  }
+  mask_ <<= iters_;
+}
+
+void SiPixelDigiMorphing::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+
+  desc.add<edm::InputTag>("src", edm::InputTag("siPixelDigis"));
+  desc.add<int32_t>("nrows", 160);
+  desc.add<int32_t>("ncols", 416);
+  desc.add<int32_t>("nrocs", 8);
+  desc.add<int32_t>("iters", 1);
+  desc.add<std::vector<int32_t>>("kernel1", {7, 7, 7});
+  desc.add<std::vector<int32_t>>("kernel2", {2, 7, 2});
+  desc.add<uint32_t>("fakeAdc", 100);
+
+  descriptions.addWithDefaultLabel(desc);
+}
+
+void SiPixelDigiMorphing::produce(edm::Event& e, const edm::EventSetup& es) {
+  auto const& inputDigi = e.get(tPixelDigi_);
+
+  auto outputDigis = std::make_unique<edm::DetSetVector<PixelDigi>>();
+
+  const int rocSize = nrows_ + 2 * iters_;
+  const int arrSize = nrocs_ * rocSize;
+
+  uint64_t imap[arrSize];
+  uint64_t map1[arrSize];
+  uint64_t map2[arrSize];
+
+  for (auto const& ds : inputDigi) {
+    auto rawId = ds.detId();
+    edm::DetSet<PixelDigi>* detDigis = nullptr;
+    detDigis = &(outputDigis->find_or_insert(rawId));
+
+    memset(imap, 0, arrSize * sizeof(uint64_t));
+    for (auto const& di : ds) {
+      int r = int(di.column()) / ncols_r_;
+      int c = int(di.column()) % ncols_r_;
+      imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c + iters_);
+      if (r > 0 && c < iters_) {
+        imap[(r - 1) * rocSize + di.row() + iters_] |= uint64_t(1) << (c + ncols_r_ + iters_);
+      } else if (++r < nrocs_ && c >= ncols_r_ - iters_) {
+        imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c - ncols_r_ + iters_);
+      }
+      (*detDigis).data.emplace_back(di.row(), di.column(), di.adc(), 0);
+    }
+
+    std::memcpy(map1, imap, arrSize * sizeof(uint64_t));
+    memset(map2, 0, arrSize * sizeof(uint64_t));
+
+    morph(map1, map2, kernel1_.data(), kDilate);
+    morph(map2, map1, kernel2_.data(), kErode);
+
+    uint64_t* i = imap + iters_;
+    uint64_t* o = map1 + iters_;
+    for (int roc = 0; roc < nrocs_; roc++, i += 2 * iters_, o += 2 * iters_) {
+      for (int row = 0; row < nrows_; row++, i++, o++) {
+        if (*o == 0)
+          continue;
+        *o >>= iters_;
+        *i >>= iters_;
+        for (int col = 0; col < ncols_r_; col++, (*i) >>= 1, (*o) >>= 1) {
+          if (*o == 0)
+            break;
+          if (((*i) & uint64_t(1)) == 1)
+            continue;
+          if (((*o) & uint64_t(1)) == 1) {
+            (*detDigis).data.emplace_back(row, roc * ncols_r_ + col, fakeAdc_, 1);
+          }
+        }
+      }
+    }
+  }
+  e.put(tPutPixelDigi_, std::move(outputDigis));
+}
+
+void SiPixelDigiMorphing::morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const {
+  uint64_t* i[ksize_];          // i(nput)
+  uint64_t* o = omap + iters_;  // o(output)
+  unsigned char valid = 0;
+  unsigned char const validMask = (1 << ksize_) - 1;
+  uint64_t m[ksize_];  // m(ask)
+
+  for (int ii = 0; ii < ksize_; ii++) {
+    i[ii] = imap + ii;
+    valid = (valid << 1) | (*i[ii] != 0);
+    m[ii] = mask_ << ii;
+  }
+
+  for (int roc = 0; roc < nrocs_; roc++, o += 2 * iters_) {
+    for (int row = 0; row < nrows_; row++, o++) {
+      if ((valid & validMask) != 0) {
+        for (int jj = 0; jj < ksize_; jj++) {
+          for (int ii = 0; ii < ksize_; ii++) {
+            uint64_t v = (*i[ii]) & (kernel[ii] << jj);  // v(ector)
+            if (op == kErode)
+              v ^= (kernel[ii] << jj);
+            uint64_t vv = v;  // vv(vector bit - shifted and contracted)
+            for (int b = 1; b < ksize_; b++)
+              vv |= (v >> b);
+            *o |= ((vv << iters_) & m[jj]);
+          }
+          if (op == kErode)
+            *o ^= m[jj];
+        }
+      }
+      for (int ii = 0; ii < ksize_; ii++)
+        i[ii]++;
+      valid = (valid << 1) | (*i[ksize_ - 1] != 0);
+    }
+    for (int ii = 0; ii < ksize_; ii++) {
+      i[ii] += 2 * iters_;
+      valid = (valid << 1) | (*i[ii] != 0);
+    }
+  }
+}
+
+DEFINE_FWK_MODULE(SiPixelDigiMorphing);

RecoLocalTracker/SiPixelDigiReProducers/python/siPixelDigisMorphed_cfi.py

@@ -0,0 +1,5 @@
+import FWCore.ParameterSet.Config as cms
+
+from RecoLocalTracker.SiPixelDigiReProducers.siPixelDigiMorphing_cfi import siPixelDigiMorphing as _siPixelDigiMorphing
+
+siPixelDigisMorphed = _siPixelDigiMorphing.clone()