SiPixelDigiReProducers/plugins/SiPixelDigiMorphing.cc

0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/stream/EDProducer.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007 #include "FWCore/PluginManager/interface/ModuleDef.h"
0008 #include "DataFormats/Common/interface/DetSetVector.h"
0009 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0010 #include "DataFormats/DetId/interface/DetId.h"
0011 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0012
0013 #include <bitset>
0014 #include <iterator>
0015
0016 class SiPixelDigiMorphing : public edm::stream::EDProducer<> {
0017 public:
0018   explicit SiPixelDigiMorphing(const edm::ParameterSet& conf);
0019   ~SiPixelDigiMorphing() override = default;
0020
0021   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0022   void produce(edm::Event& e, const edm::EventSetup& c) override;
0023
0024 private:
0025   edm::EDGetTokenT<edm::DetSetVector<PixelDigi>> tPixelDigi_;
0026   edm::EDPutTokenT<edm::DetSetVector<PixelDigi>> tPutPixelDigi_;
0027
0028   const int32_t nrows_;
0029   const int32_t ncols_;
0030   const int32_t nrocs_;  // in Phase 1, this is ROCs, but could be any subset of a pixel row
0031   const int32_t iters_;
0032   const uint32_t fakeAdc_;
0033
0034   int32_t ncols_r_;  // number of columns per ROC
0035   int32_t ksize_;    // kernel size
0036
0037   std::vector<uint64_t> kernel1_;
0038   std::vector<uint64_t> kernel2_;
0039   uint64_t mask_;
0040
0041   enum MorphOption { kDilate, kErode };
0042
0043   void morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const;
0044 };
0045
0046 SiPixelDigiMorphing::SiPixelDigiMorphing(edm::ParameterSet const& conf)
0047     : tPixelDigi_(consumes(conf.getParameter<edm::InputTag>("src"))),
0048       tPutPixelDigi_(produces<edm::DetSetVector<PixelDigi>>()),
0049       nrows_(conf.getParameter<int32_t>("nrows")),
0050       ncols_(conf.getParameter<int32_t>("ncols")),
0051       nrocs_(conf.getParameter<int32_t>("nrocs")),
0052       iters_(conf.getParameter<int32_t>("iters")),
0053       fakeAdc_(conf.getParameter<uint32_t>("fakeAdc")) {
0054   if (ncols_ % nrocs_ == 0) {
0055     ncols_r_ = ncols_ / nrocs_;
0056   } else {
0057     throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
0058                                           << " number of columns not divisible with"
0059                                           << " number of ROCs\n";
0060   }
0061
0062   if (ncols_r_ + 2 * iters_ <= int(sizeof(uint64_t) * 8)) {
0063     ksize_ = 2 * iters_ + 1;
0064   } else {
0065     throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
0066                                           << " too many columns per ROC"
0067                                           << " or too many iterations set\n"
0068                                           << " Ncol/Nrocs+2*iters should not be"
0069                                           << " more than " << sizeof(uint64_t) * 8 << "\n";
0070   }
0071
0072   std::vector<int32_t> k1(conf.getParameter<std::vector<int32_t>>("kernel1"));
0073   std::vector<int32_t> k2(conf.getParameter<std::vector<int32_t>>("kernel2"));
0074
0075   kernel1_.resize(ksize_, 0);
0076   kernel2_.resize(ksize_, 0);
0077   mask_ = 0;
0078   int w = (ncols_r_ + 2 * iters_) / ksize_ + ((ncols_r_ + 2 * iters_) % ksize_ != 0);
0079   for (int j = 0; j < w; j++) {
0080     for (int ii = 0; ii < ksize_; ii++) {
0081       kernel1_[ii] <<= ksize_;
0082       kernel1_[ii] |= k1[ii];
0083       kernel2_[ii] <<= ksize_;
0084       kernel2_[ii] |= k2[ii];
0085     }
0086     mask_ <<= ksize_;
0087     mask_ |= 1;
0088   }
0089   mask_ <<= iters_;
0090 }
0091
0092 void SiPixelDigiMorphing::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0093   edm::ParameterSetDescription desc;
0094
0095   desc.add<edm::InputTag>("src", edm::InputTag("siPixelDigis"));
0096   desc.add<int32_t>("nrows", 160);
0097   desc.add<int32_t>("ncols", 416);
0098   desc.add<int32_t>("nrocs", 8);
0099   desc.add<int32_t>("iters", 1);
0100   desc.add<std::vector<int32_t>>("kernel1", {7, 7, 7});
0101   desc.add<std::vector<int32_t>>("kernel2", {2, 7, 2});
0102   desc.add<uint32_t>("fakeAdc", 100);
0103
0104   descriptions.addWithDefaultLabel(desc);
0105 }
0106
0107 void SiPixelDigiMorphing::produce(edm::Event& e, const edm::EventSetup& es) {
0108   auto const& inputDigi = e.get(tPixelDigi_);
0109
0110   auto outputDigis = std::make_unique<edm::DetSetVector<PixelDigi>>();
0111
0112   const int rocSize = nrows_ + 2 * iters_;
0113   const int arrSize = nrocs_ * rocSize;
0114
0115   uint64_t imap[arrSize];
0116   uint64_t map1[arrSize];
0117   uint64_t map2[arrSize];
0118
0119   for (auto const& ds : inputDigi) {
0120     auto rawId = ds.detId();
0121     edm::DetSet<PixelDigi>* detDigis = nullptr;
0122     detDigis = &(outputDigis->find_or_insert(rawId));
0123
0124     memset(imap, 0, arrSize * sizeof(uint64_t));
0125     for (auto const& di : ds) {
0126       int r = int(di.column()) / ncols_r_;
0127       int c = int(di.column()) % ncols_r_;
0128       imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c + iters_);
0129       if (r > 0 && c < iters_) {
0130         imap[(r - 1) * rocSize + di.row() + iters_] |= uint64_t(1) << (c + ncols_r_ + iters_);
0131       } else if (++r < nrocs_ && c >= ncols_r_ - iters_) {
0132         imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c - ncols_r_ + iters_);
0133       }
0134       (*detDigis).data.emplace_back(di.row(), di.column(), di.adc(), 0);
0135     }
0136
0137     std::memcpy(map1, imap, arrSize * sizeof(uint64_t));
0138     memset(map2, 0, arrSize * sizeof(uint64_t));
0139
0140     morph(map1, map2, kernel1_.data(), kDilate);
0141     morph(map2, map1, kernel2_.data(), kErode);
0142
0143     uint64_t* i = imap + iters_;
0144     uint64_t* o = map1 + iters_;
0145     for (int roc = 0; roc < nrocs_; roc++, i += 2 * iters_, o += 2 * iters_) {
0146       for (int row = 0; row < nrows_; row++, i++, o++) {
0147         if (*o == 0)
0148           continue;
0149         *o >>= iters_;
0150         *i >>= iters_;
0151         for (int col = 0; col < ncols_r_; col++, (*i) >>= 1, (*o) >>= 1) {
0152           if (*o == 0)
0153             break;
0154           if (((*i) & uint64_t(1)) == 1)
0155             continue;
0156           if (((*o) & uint64_t(1)) == 1) {
0157             (*detDigis).data.emplace_back(row, roc * ncols_r_ + col, fakeAdc_, 1);
0158           }
0159         }
0160       }
0161     }
0162   }
0163   e.put(tPutPixelDigi_, std::move(outputDigis));
0164 }
0165
0166 void SiPixelDigiMorphing::morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const {
0167   uint64_t* i[ksize_];          // i(nput)
0168   uint64_t* o = omap + iters_;  // o(output)
0169   unsigned char valid = 0;
0170   unsigned char const validMask = (1 << ksize_) - 1;
0171   uint64_t m[ksize_];  // m(ask)
0172
0173   for (int ii = 0; ii < ksize_; ii++) {
0174     i[ii] = imap + ii;
0175     valid = (valid << 1) | (*i[ii] != 0);
0176     m[ii] = mask_ << ii;
0177   }
0178
0179   for (int roc = 0; roc < nrocs_; roc++, o += 2 * iters_) {
0180     for (int row = 0; row < nrows_; row++, o++) {
0181       if ((valid & validMask) != 0) {
0182         for (int jj = 0; jj < ksize_; jj++) {
0183           for (int ii = 0; ii < ksize_; ii++) {
0184             uint64_t v = (*i[ii]) & (kernel[ii] << jj);  // v(ector)
0185             if (op == kErode)
0186               v ^= (kernel[ii] << jj);
0187             uint64_t vv = v;  // vv(vector bit - shifted and contracted)
0188             for (int b = 1; b < ksize_; b++)
0189               vv |= (v >> b);
0190             *o |= ((vv << iters_) & m[jj]);
0191           }
0192           if (op == kErode)
0193             *o ^= m[jj];
0194         }
0195       }
0196       for (int ii = 0; ii < ksize_; ii++)
0197         i[ii]++;
0198       valid = (valid << 1) | (*i[ksize_ - 1] != 0);
0199     }
0200     for (int ii = 0; ii < ksize_; ii++) {
0201       i[ii] += 2 * iters_;
0202       valid = (valid << 1) | (*i[ii] != 0);
0203     }
0204   }
0205 }
0206
0207 DEFINE_FWK_MODULE(SiPixelDigiMorphing);