SiPixelClusterizer/plugins/SiPixelDigisClustersFromSoAAlpaka.cc

0001 #include <memory>
0002
0003 #include "DataFormats/Common/interface/DetSetVector.h"
0004 #include "DataFormats/DetId/interface/DetId.h"
0005 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0006 #include "DataFormats/SiPixelClusterSoA/interface/ClusteringConstants.h"
0007 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0008 #include "DataFormats/SiPixelDigiSoA/interface/SiPixelDigisHost.h"
0009 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/global/EDProducer.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0017 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0018 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0019 #include "RecoLocalTracker/SiPixelClusterizer/interface/SiPixelClusterThresholds.h"
0020
0021 // local include(s)
0022 #include "PixelClusterizerBase.h"
0023
0024 //#define EDM_ML_DEBUG
0025 //#define GPU_DEBUG
0026
0027 template <typename TrackerTraits>
0028 class SiPixelDigisClustersFromSoAAlpaka : public edm::global::EDProducer<> {
0029 public:
0030   explicit SiPixelDigisClustersFromSoAAlpaka(const edm::ParameterSet& iConfig);
0031   ~SiPixelDigisClustersFromSoAAlpaka() override = default;
0032
0033   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0034
0035 private:
0036   void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0037
0038   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> const topoToken_;
0039   edm::EDGetTokenT<SiPixelDigisHost> const digisHostToken_;
0040   const SiPixelClusterThresholds clusterThresholds_;  // Cluster threshold in electrons
0041   const bool produceDigis_;
0042   const bool storeDigis_;
0043
0044   edm::EDPutTokenT<edm::DetSetVector<PixelDigi>> digisPutToken_;
0045   edm::EDPutTokenT<SiPixelClusterCollectionNew> clustersPutToken_;
0046 };
0047
0048 template <typename TrackerTraits>
0049 SiPixelDigisClustersFromSoAAlpaka<TrackerTraits>::SiPixelDigisClustersFromSoAAlpaka(const edm::ParameterSet& iConfig)
0050     : topoToken_(esConsumes()),
0051       digisHostToken_(consumes(iConfig.getParameter<edm::InputTag>("src"))),
0052       clusterThresholds_(iConfig.getParameter<int>("clusterThreshold_layer1"),
0053                          iConfig.getParameter<int>("clusterThreshold_otherLayers")),
0054       produceDigis_(iConfig.getParameter<bool>("produceDigis")),
0055       storeDigis_(produceDigis_ && iConfig.getParameter<bool>("storeDigis")),
0056       clustersPutToken_(produces()) {
0057   if (produceDigis_)
0058     digisPutToken_ = produces<edm::DetSetVector<PixelDigi>>();
0059 }
0060
0061 template <typename TrackerTraits>
0062 void SiPixelDigisClustersFromSoAAlpaka<TrackerTraits>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0063   edm::ParameterSetDescription desc;
0064   desc.add<edm::InputTag>("src", edm::InputTag("siPixelDigisSoA"));
0065   desc.add<int>("clusterThreshold_layer1", pixelClustering::clusterThresholdLayerOne);
0066   desc.add<int>("clusterThreshold_otherLayers", pixelClustering::clusterThresholdOtherLayers);
0067   desc.add<bool>("produceDigis", true);
0068   desc.add<bool>("storeDigis", true);
0069
0070   descriptions.addWithDefaultLabel(desc);
0071 }
0072
0073 template <typename TrackerTraits>
0074 void SiPixelDigisClustersFromSoAAlpaka<TrackerTraits>::produce(edm::StreamID,
0075                                                                edm::Event& iEvent,
0076                                                                const edm::EventSetup& iSetup) const {
0077   const auto& digisHost = iEvent.get(digisHostToken_);
0078   const auto& digisView = digisHost.const_view();
0079   const uint32_t nDigis = digisHost.nDigis();
0080
0081   const auto& ttopo = iSetup.getData(topoToken_);
0082   constexpr auto maxModules = TrackerTraits::numberOfModules;
0083
0084   std::unique_ptr<edm::DetSetVector<PixelDigi>> outputDigis;
0085   if (produceDigis_)
0086     outputDigis = std::make_unique<edm::DetSetVector<PixelDigi>>();
0087   if (storeDigis_)
0088     outputDigis->reserve(maxModules);
0089   auto outputClusters = std::make_unique<SiPixelClusterCollectionNew>();
0090   outputClusters->reserve(maxModules, nDigis / 2);
0091
0092   edm::DetSet<PixelDigi>* detDigis = nullptr;
0093   uint32_t detId = 0;
0094
0095   for (uint32_t i = 0; i < nDigis; i++) {
0096     // check for uninitialized digis
0097     // this is set in RawToDigi_kernel in SiPixelRawToClusterGPUKernel.cu
0098     if (digisView[i].rawIdArr() == 0)
0099       continue;
0100
0101     // check for noisy/dead pixels (electrons set to 0)
0102     if (digisView[i].adc() == 0)
0103       continue;
0104
0105     detId = digisView[i].rawIdArr();
0106     if (storeDigis_) {
0107       detDigis = &outputDigis->find_or_insert(detId);
0108
0109       if ((*detDigis).empty())
0110         (*detDigis).data.reserve(64);  // avoid the first relocations
0111     }
0112
0113     break;
0114   }
0115
0116   int32_t nclus = -1;
0117   PixelClusterizerBase::AccretionCluster aclusters[TrackerTraits::maxNumClustersPerModules];
0118 #ifdef EDM_ML_DEBUG
0119   auto totClustersFilled = 0;
0120 #endif
0121
0122   auto fillClusters = [&](uint32_t detId) {
0123     if (nclus < 0)
0124       return;  // this in reality should never happen
0125     edmNew::DetSetVector<SiPixelCluster>::FastFiller spc(*outputClusters, detId);
0126     auto layer = (DetId(detId).subdetId() == 1) ? ttopo.pxbLayer(detId) : 0;
0127     auto clusterThreshold = clusterThresholds_.getThresholdForLayerOnCondition(layer == 1);
0128     for (int32_t ic = 0; ic < nclus + 1; ++ic) {
0129       auto const& acluster = aclusters[ic];
0130       // in any case we cannot  go out of sync with gpu...
0131       if (acluster.charge < clusterThreshold)
0132         edm::LogWarning("SiPixelDigisClustersFromSoAAlpaka")
0133             << "cluster below charge Threshold "
0134             << "Layer/DetId/clusId " << layer << '/' << detId << '/' << ic << " size/charge " << acluster.isize << '/'
0135             << acluster.charge << "\n";
0136       // sort by row (x)
0137       spc.emplace_back(acluster.isize, acluster.adc, acluster.x, acluster.y, acluster.xmin, acluster.ymin, ic);
0138       aclusters[ic].clear();
0139 #ifdef EDM_ML_DEBUG
0140       ++totClustersFilled;
0141       const auto& cluster{spc.back()};
0142       // LogDebug("SiPixelDigisClustersFromSoAAlpaka")
0143       std::cout << "putting in this cluster " << ic << " " << cluster.charge() << " " << cluster.pixelADC().size()
0144                 << "\n";
0145 #endif
0146       std::push_heap(spc.begin(), spc.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0147         return cl1.minPixelRow() < cl2.minPixelRow();
0148       });
0149     }
0150     nclus = -1;
0151     // sort by row (x)
0152     std::sort_heap(spc.begin(), spc.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0153       return cl1.minPixelRow() < cl2.minPixelRow();
0154     });
0155     if (spc.empty())
0156       spc.abort();
0157   };
0158
0159 #ifdef GPU_DEBUG
0160   std::cout << "Dumping all digis. nDigis = " << nDigis << std::endl;
0161 #endif
0162   for (uint32_t i = 0; i < nDigis; i++) {
0163     // check for uninitialized digis
0164     if (digisView[i].rawIdArr() == 0)
0165       continue;
0166     // check for noisy/dead pixels (electrons set to 0)
0167     if (digisView[i].adc() == 0)
0168       continue;
0169     // not in cluster; TODO add an assert for the size
0170     if (digisView[i].clus() == pixelClustering::invalidClusterId) {
0171       continue;
0172     }
0173     // unexpected invalid value
0174     if (digisView[i].clus() < pixelClustering::invalidClusterId) {
0175       edm::LogError("SiPixelDigisClustersFromSoAAlpaka")
0176           << "Skipping pixel digi with unexpected invalid cluster id " << digisView[i].clus();
0177       continue;
0178     }
0179     // from clusters killed by charge cut
0180     if (digisView[i].clus() == pixelClustering::invalidModuleId)
0181       continue;
0182
0183 #ifdef EDM_ML_DEBUG
0184     assert(digisView[i].rawIdArr() > 109999);
0185 #endif
0186     if (detId != digisView[i].rawIdArr()) {
0187 #ifdef GPU_DEBUG
0188       std::cout << ">> Closed module --" << detId << "; nclus = " << nclus << std::endl;
0189 #endif
0190       // new module
0191       fillClusters(detId);
0192 #ifdef EDM_ML_DEBUG
0193       assert(nclus == -1);
0194 #endif
0195       detId = digisView[i].rawIdArr();
0196       if (storeDigis_) {
0197         detDigis = &outputDigis->find_or_insert(detId);
0198         if ((*detDigis).empty())
0199           (*detDigis).data.reserve(64);  // avoid the first relocations
0200         else {
0201           edm::LogWarning("SiPixelDigisClustersFromSoAAlpaka")
0202               << "Problem det present twice in input! " << (*detDigis).detId();
0203         }
0204       }
0205     }
0206     PixelDigi dig{digisView[i].pdigi()};
0207 #ifdef GPU_DEBUG
0208     std::cout << i << ";" << digisView[i].rawIdArr() << ";" << digisView[i].clus() << ";" << digisView[i].pdigi() << ";"
0209               << digisView[i].adc() << ";" << dig.row() << ";" << dig.column() << std::endl;
0210 #endif
0211
0212     if (storeDigis_)
0213       (*detDigis).data.emplace_back(dig);
0214       // fill clusters
0215 #ifdef EDM_ML_DEBUG
0216     assert(digisView[i].clus() >= 0);
0217     assert(digisView[i].clus() < static_cast<int>(TrackerTraits::maxNumClustersPerModules));
0218 #endif
0219     nclus = std::max(digisView[i].clus(), nclus);
0220     auto row = dig.row();
0221     auto col = dig.column();
0222     SiPixelCluster::PixelPos pix(row, col);
0223     aclusters[digisView[i].clus()].add(pix, digisView[i].adc());
0224   }
0225
0226   // fill final clusters
0227   if (detId > 0)
0228     fillClusters(detId);
0229
0230 #ifdef EDM_ML_DEBUG
0231   LogDebug("SiPixelDigisClustersFromSoAAlpaka") << "filled " << totClustersFilled << " clusters";
0232 #endif
0233
0234   if (produceDigis_)
0235     iEvent.put(digisPutToken_, std::move(outputDigis));
0236
0237   iEvent.put(clustersPutToken_, std::move(outputClusters));
0238 }
0239
0240 #include "FWCore/Framework/interface/MakerMacros.h"
0241
0242 using SiPixelDigisClustersFromSoAAlpakaPhase1 = SiPixelDigisClustersFromSoAAlpaka<pixelTopology::Phase1>;
0243 DEFINE_FWK_MODULE(SiPixelDigisClustersFromSoAAlpakaPhase1);
0244
0245 using SiPixelDigisClustersFromSoAAlpakaPhase2 = SiPixelDigisClustersFromSoAAlpaka<pixelTopology::Phase2>;
0246 DEFINE_FWK_MODULE(SiPixelDigisClustersFromSoAAlpakaPhase2);
0247
0248 using SiPixelDigisClustersFromSoAAlpakaHIonPhase1 = SiPixelDigisClustersFromSoAAlpaka<pixelTopology::HIonPhase1>;
0249 DEFINE_FWK_MODULE(SiPixelDigisClustersFromSoAAlpakaHIonPhase1);