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 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0013 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0014 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0015
0016 #include <boost/algorithm/string.hpp>
0017 #include <bitset>
0018 #include <iterator>
0019
0020 typedef std::pair<uint32_t, uint32_t> range;
0021 typedef std::vector<range> region;
0022
0023 class SiPixelDigiMorphing : public edm::stream::EDProducer<> {
0024 public:
0025   explicit SiPixelDigiMorphing(const edm::ParameterSet& conf);
0026   ~SiPixelDigiMorphing() override = default;
0027
0028   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0029   void produce(edm::Event& e, const edm::EventSetup& c) override;
0030
0031 private:
0032   edm::EDGetTokenT<edm::DetSetVector<PixelDigi>> tPixelDigi_;
0033   edm::EDPutTokenT<edm::DetSetVector<PixelDigi>> tPutPixelDigi_;
0034   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopologyToken_;
0035
0036   const std::vector<region> theBarrelRegions_;
0037   const std::vector<region> theEndcapRegions_;
0038
0039   const int32_t nrows_;
0040   const int32_t ncols_;
0041   const int32_t nrocs_;  // in Phase 1, this is ROCs, but could be any subset of a pixel row
0042   const int32_t iters_;
0043   const uint32_t fakeAdc_;
0044
0045   int32_t ncols_r_;  // number of columns per ROC
0046   int32_t ksize_;    // kernel size
0047
0048   std::vector<uint64_t> kernel1_;
0049   std::vector<uint64_t> kernel2_;
0050   uint64_t mask_;
0051
0052   enum MorphOption { kDilate, kErode };
0053
0054   void morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const;
0055   std::vector<region> parseRegions(const std::vector<std::string>& regionStrings, size_t size);
0056   bool skipDetId(const TrackerTopology* tTopo,
0057                  const DetId& detId,
0058                  const std::vector<region>& theBarrelRegions,
0059                  const std::vector<region>& theEndcapRegions) const;
0060 };
0061
0062 SiPixelDigiMorphing::SiPixelDigiMorphing(edm::ParameterSet const& conf)
0063     : tPixelDigi_(consumes(conf.getParameter<edm::InputTag>("src"))),
0064       tPutPixelDigi_(produces<edm::DetSetVector<PixelDigi>>()),
0065       trackerTopologyToken_(esConsumes()),
0066       theBarrelRegions_(parseRegions(conf.getParameter<std::vector<std::string>>("barrelRegions"), 3)),
0067       theEndcapRegions_(parseRegions(conf.getParameter<std::vector<std::string>>("endcapRegions"), 4)),
0068       nrows_(conf.getParameter<int32_t>("nrows")),
0069       ncols_(conf.getParameter<int32_t>("ncols")),
0070       nrocs_(conf.getParameter<int32_t>("nrocs")),
0071       iters_(conf.getParameter<int32_t>("iters")),
0072       fakeAdc_(conf.getParameter<uint32_t>("fakeAdc")) {
0073   if (ncols_ % nrocs_ == 0) {
0074     ncols_r_ = ncols_ / nrocs_;
0075   } else {
0076     throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
0077                                           << " number of columns not divisible with"
0078                                           << " number of ROCs\n";
0079   }
0080
0081   if (ncols_r_ + 2 * iters_ <= int(sizeof(uint64_t) * 8)) {
0082     ksize_ = 2 * iters_ + 1;
0083   } else {
0084     throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
0085                                           << " too many columns per ROC"
0086                                           << " or too many iterations set\n"
0087                                           << " Ncol/Nrocs+2*iters should not be"
0088                                           << " more than " << sizeof(uint64_t) * 8 << "\n";
0089   }
0090
0091   std::vector<int32_t> k1(conf.getParameter<std::vector<int32_t>>("kernel1"));
0092   std::vector<int32_t> k2(conf.getParameter<std::vector<int32_t>>("kernel2"));
0093
0094   kernel1_.resize(ksize_, 0);
0095   kernel2_.resize(ksize_, 0);
0096   mask_ = 0;
0097   int w = (ncols_r_ + 2 * iters_) / ksize_ + ((ncols_r_ + 2 * iters_) % ksize_ != 0);
0098   for (int j = 0; j < w; j++) {
0099     for (int ii = 0; ii < ksize_; ii++) {
0100       kernel1_[ii] <<= ksize_;
0101       kernel1_[ii] |= k1[ii];
0102       kernel2_[ii] <<= ksize_;
0103       kernel2_[ii] |= k2[ii];
0104     }
0105     mask_ <<= ksize_;
0106     mask_ |= 1;
0107   }
0108   mask_ <<= iters_;
0109 }
0110
0111 void SiPixelDigiMorphing::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0112   edm::ParameterSetDescription desc;
0113
0114   desc.add<edm::InputTag>("src", edm::InputTag("siPixelDigis"));
0115   // LAYER,LADDER,MODULE (coordinates can also be specified as a range FIRST-LAST where appropriate)
0116   desc.add<std::vector<std::string>>("barrelRegions", {"1,1-12,1-2", "1,1-12,7-8", "2,1-28,1", "1,1-28,8"});
0117   // DISK,BLADE,SIDE,PANEL (coordinates can also be specified as a range FIRST-LAST where appropriate)
0118   desc.add<std::vector<std::string>>("endcapRegions", {});
0119   desc.add<int32_t>("nrows", 160);
0120   desc.add<int32_t>("ncols", 416);
0121   desc.add<int32_t>("nrocs", 8);
0122   desc.add<int32_t>("iters", 1);
0123   desc.add<std::vector<int32_t>>("kernel1", {7, 7, 7});
0124   desc.add<std::vector<int32_t>>("kernel2", {2, 7, 2});
0125   desc.add<uint32_t>("fakeAdc", 100);
0126
0127   descriptions.addWithDefaultLabel(desc);
0128 }
0129
0130 void SiPixelDigiMorphing::produce(edm::Event& e, const edm::EventSetup& es) {
0131   auto const& inputDigi = e.get(tPixelDigi_);
0132
0133   auto outputDigis = std::make_unique<edm::DetSetVector<PixelDigi>>();
0134
0135   // use the TrackerTopology for the layer/disk number, etc.
0136   const TrackerTopology* tTopo = &es.getData(trackerTopologyToken_);
0137
0138   const int rocSize = nrows_ + 2 * iters_;
0139   const int arrSize = nrocs_ * rocSize;
0140
0141   uint64_t imap[arrSize];
0142   uint64_t map1[arrSize];
0143   uint64_t map2[arrSize];
0144
0145   for (auto const& ds : inputDigi) {
0146     auto rawId = ds.detId();
0147
0148     const DetId detId(rawId);
0149
0150     // skip DetIds for which digi morphing has not been requested
0151     if (skipDetId(tTopo, detId, theBarrelRegions_, theEndcapRegions_)) {
0152       outputDigis->insert(ds);
0153       continue;
0154     }
0155
0156     edm::DetSet<PixelDigi>* detDigis = nullptr;
0157     detDigis = &(outputDigis->find_or_insert(rawId));
0158
0159     memset(imap, 0, arrSize * sizeof(uint64_t));
0160     for (auto const& di : ds) {
0161       int r = int(di.column()) / ncols_r_;
0162       int c = int(di.column()) % ncols_r_;
0163       imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c + iters_);
0164       if (r > 0 && c < iters_) {
0165         imap[(r - 1) * rocSize + di.row() + iters_] |= uint64_t(1) << (c + ncols_r_ + iters_);
0166       } else if (++r < nrocs_ && c >= ncols_r_ - iters_) {
0167         imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c - ncols_r_ + iters_);
0168       }
0169       (*detDigis).data.emplace_back(di.row(), di.column(), di.adc(), 0);
0170     }
0171
0172     std::memcpy(map1, imap, arrSize * sizeof(uint64_t));
0173     memset(map2, 0, arrSize * sizeof(uint64_t));
0174
0175     morph(map1, map2, kernel1_.data(), kDilate);
0176     morph(map2, map1, kernel2_.data(), kErode);
0177
0178     uint64_t* i = imap + iters_;
0179     uint64_t* o = map1 + iters_;
0180     for (int roc = 0; roc < nrocs_; roc++, i += 2 * iters_, o += 2 * iters_) {
0181       for (int row = 0; row < nrows_; row++, i++, o++) {
0182         if (*o == 0)
0183           continue;
0184         *o >>= iters_;
0185         *i >>= iters_;
0186         for (int col = 0; col < ncols_r_; col++, (*i) >>= 1, (*o) >>= 1) {
0187           if (*o == 0)
0188             break;
0189           if (((*i) & uint64_t(1)) == 1)
0190             continue;
0191           if (((*o) & uint64_t(1)) == 1) {
0192             (*detDigis).data.emplace_back(row, roc * ncols_r_ + col, fakeAdc_, 1);
0193           }
0194         }
0195       }
0196     }
0197   }
0198   e.put(tPutPixelDigi_, std::move(outputDigis));
0199 }
0200
0201 void SiPixelDigiMorphing::morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const {
0202   uint64_t* i[ksize_];          // i(nput)
0203   uint64_t* o = omap + iters_;  // o(output)
0204   unsigned char valid = 0;
0205   unsigned char const validMask = (1 << ksize_) - 1;
0206   uint64_t m[ksize_];  // m(ask)
0207
0208   for (int ii = 0; ii < ksize_; ii++) {
0209     i[ii] = imap + ii;
0210     valid = (valid << 1) | (*i[ii] != 0);
0211     m[ii] = mask_ << ii;
0212   }
0213
0214   for (int roc = 0; roc < nrocs_; roc++, o += 2 * iters_) {
0215     for (int row = 0; row < nrows_; row++, o++) {
0216       if ((valid & validMask) != 0) {
0217         for (int jj = 0; jj < ksize_; jj++) {
0218           for (int ii = 0; ii < ksize_; ii++) {
0219             uint64_t v = (*i[ii]) & (kernel[ii] << jj);  // v(ector)
0220             if (op == kErode)
0221               v ^= (kernel[ii] << jj);
0222             uint64_t vv = v;  // vv(vector bit - shifted and contracted)
0223             for (int b = 1; b < ksize_; b++)
0224               vv |= (v >> b);
0225             *o |= ((vv << iters_) & m[jj]);
0226           }
0227           if (op == kErode)
0228             *o ^= m[jj];
0229         }
0230       }
0231       for (int ii = 0; ii < ksize_; ii++)
0232         i[ii]++;
0233       valid = (valid << 1) | (*i[ksize_ - 1] != 0);
0234     }
0235     for (int ii = 0; ii < ksize_; ii++) {
0236       i[ii] += 2 * iters_;
0237       valid = (valid << 1) | (*i[ii] != 0);
0238     }
0239   }
0240 }
0241
0242 std::vector<region> SiPixelDigiMorphing::parseRegions(const std::vector<std::string>& regionStrings, size_t size) {
0243   std::vector<region> regions;
0244
0245   for (auto const& str : regionStrings) {
0246     region reg;
0247
0248     std::vector<std::string> ranges;
0249     boost::split(ranges, str, boost::is_any_of(","));
0250
0251     if (ranges.size() != size) {
0252       throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
0253                                             << " invalid number of coordinates provided in " << str << " (" << size
0254                                             << " expected, " << ranges.size() << " provided)\n";
0255     }
0256
0257     for (auto const& r : ranges) {
0258       std::vector<std::string> limits;
0259       boost::split(limits, r, boost::is_any_of("-"));
0260
0261       try {
0262         // if range specified
0263         if (limits.size() > 1) {
0264           reg.push_back(std::make_pair(std::stoi(limits.at(0)), std::stoi(limits.at(1))));
0265           // otherwise store single value as a range
0266         } else {
0267           reg.push_back(std::make_pair(std::stoi(limits.at(0)), std::stoi(limits.at(0))));
0268         }
0269       } catch (...) {
0270         throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
0271                                               << " invalid coordinate value provided in " << str << "\n";
0272       }
0273     }
0274     regions.push_back(reg);
0275   }
0276
0277   return regions;
0278 }
0279
0280 // apply regional digi morphing logic
0281 bool SiPixelDigiMorphing::skipDetId(const TrackerTopology* tTopo,
0282                                     const DetId& detId,
0283                                     const std::vector<region>& theBarrelRegions,
0284                                     const std::vector<region>& theEndcapRegions) const {
0285   // barrel
0286   if (detId.subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel)) {
0287     // no barrel region specified
0288     if (theBarrelRegions.empty()) {
0289       return true;
0290     } else {
0291       uint32_t layer = tTopo->pxbLayer(detId.rawId());
0292       uint32_t ladder = tTopo->pxbLadder(detId.rawId());
0293       uint32_t module = tTopo->pxbModule(detId.rawId());
0294
0295       bool inRegion = false;
0296
0297       for (auto const& reg : theBarrelRegions) {
0298         if ((layer >= reg.at(0).first && layer <= reg.at(0).second) &&
0299             (ladder >= reg.at(1).first && ladder <= reg.at(1).second) &&
0300             (module >= reg.at(2).first && module <= reg.at(2).second)) {
0301           inRegion = true;
0302           break;
0303         }
0304       }
0305
0306       return !inRegion;
0307     }
0308     // endcap
0309   } else {
0310     // no endcap region specified
0311     if (theEndcapRegions.empty()) {
0312       return true;
0313     } else {
0314       uint32_t disk = tTopo->pxfDisk(detId.rawId());
0315       uint32_t blade = tTopo->pxfBlade(detId.rawId());
0316       uint32_t side = tTopo->pxfSide(detId.rawId());
0317       uint32_t panel = tTopo->pxfPanel(detId.rawId());
0318
0319       bool inRegion = false;
0320
0321       for (auto const& reg : theEndcapRegions) {
0322         if ((disk >= reg.at(0).first && disk <= reg.at(0).second) &&
0323             (blade >= reg.at(1).first && blade <= reg.at(1).second) &&
0324             (side >= reg.at(2).first && side <= reg.at(2).second) &&
0325             (panel >= reg.at(3).first && panel <= reg.at(3).second)) {
0326           inRegion = true;
0327           break;
0328         }
0329       }
0330
0331       return !inRegion;
0332     }
0333   }
0334 }
0335
0336 DEFINE_FWK_MODULE(SiPixelDigiMorphing);