SiPixelClusterizer/plugins/SiPixelDigisClustersFromSoA.cc

0001 #include "CUDADataFormats/SiPixelCluster/interface/gpuClusteringConstants.h"
0002 #include "DataFormats/Common/interface/DetSetVector.h"
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "DataFormats/DetId/interface/DetId.h"
0005 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0006 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0007 #include "DataFormats/SiPixelDigi/interface/SiPixelDigisSoA.h"
0008 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/MakerMacros.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/Records/interface/TrackerTopologyRcd.h"
0018 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0019
0020 // local include(s)
0021 #include "PixelClusterizerBase.h"
0022 #include "SiPixelClusterThresholds.h"
0023
0024 template <typename TrackerTraits>
0025 class SiPixelDigisClustersFromSoAT : public edm::global::EDProducer<> {
0026 public:
0027   explicit SiPixelDigisClustersFromSoAT(const edm::ParameterSet& iConfig);
0028   ~SiPixelDigisClustersFromSoAT() override = default;
0029
0030   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0031
0032 private:
0033   void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0034
0035   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0036
0037   edm::EDGetTokenT<SiPixelDigisSoA> digiGetToken_;
0038
0039   edm::EDPutTokenT<edm::DetSetVector<PixelDigi>> digiPutToken_;
0040   edm::EDPutTokenT<SiPixelClusterCollectionNew> clusterPutToken_;
0041
0042   const SiPixelClusterThresholds clusterThresholds_;  // Cluster threshold in electrons
0043
0044   const bool produceDigis_;
0045   const bool storeDigis_;
0046 };
0047
0048 template <typename TrackerTraits>
0049 SiPixelDigisClustersFromSoAT<TrackerTraits>::SiPixelDigisClustersFromSoAT(const edm::ParameterSet& iConfig)
0050     : topoToken_(esConsumes()),
0051       digiGetToken_(consumes<SiPixelDigisSoA>(iConfig.getParameter<edm::InputTag>("src"))),
0052       clusterPutToken_(produces<SiPixelClusterCollectionNew>()),
0053       clusterThresholds_(iConfig.getParameter<int>("clusterThreshold_layer1"),
0054                          iConfig.getParameter<int>("clusterThreshold_otherLayers")),
0055       produceDigis_(iConfig.getParameter<bool>("produceDigis")),
0056       storeDigis_(iConfig.getParameter<bool>("produceDigis") && iConfig.getParameter<bool>("storeDigis")) {
0057   if (produceDigis_)
0058     digiPutToken_ = produces<edm::DetSetVector<PixelDigi>>();
0059 }
0060
0061 template <typename TrackerTraits>
0062 void SiPixelDigisClustersFromSoAT<TrackerTraits>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0063   edm::ParameterSetDescription desc;
0064   desc.add<edm::InputTag>("src", edm::InputTag("siPixelDigisSoA"));
0065   desc.add<int>("clusterThreshold_layer1", gpuClustering::clusterThresholdLayerOne);
0066   desc.add<int>("clusterThreshold_otherLayers", gpuClustering::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 SiPixelDigisClustersFromSoAT<TrackerTraits>::produce(edm::StreamID,
0075                                                           edm::Event& iEvent,
0076                                                           const edm::EventSetup& iSetup) const {
0077   const auto& digis = iEvent.get(digiGetToken_);
0078   const uint32_t nDigis = digis.size();
0079   const auto& ttopo = iSetup.getData(topoToken_);
0080   constexpr auto maxModules = TrackerTraits::numberOfModules;
0081
0082   std::unique_ptr<edm::DetSetVector<PixelDigi>> collection;
0083   if (produceDigis_)
0084     collection = std::make_unique<edm::DetSetVector<PixelDigi>>();
0085   if (storeDigis_)
0086     collection->reserve(maxModules);
0087   auto outputClusters = std::make_unique<SiPixelClusterCollectionNew>();
0088   outputClusters->reserve(maxModules, nDigis / 2);
0089
0090   edm::DetSet<PixelDigi>* detDigis = nullptr;
0091   uint32_t detId = 0;
0092   for (uint32_t i = 0; i < nDigis; i++) {
0093     // check for uninitialized digis
0094     // this is set in RawToDigi_kernel in SiPixelRawToClusterGPUKernel.cu
0095     if (digis.rawIdArr(i) == 0)
0096       continue;
0097     // check for noisy/dead pixels (electrons set to 0)
0098     if (digis.adc(i) == 0)
0099       continue;
0100
0101     detId = digis.rawIdArr(i);
0102     if (storeDigis_) {
0103       detDigis = &collection->find_or_insert(detId);
0104       if ((*detDigis).empty())
0105         (*detDigis).data.reserve(64);  // avoid the first relocations
0106     }
0107     break;
0108   }
0109
0110   int32_t nclus = -1;
0111   PixelClusterizerBase::AccretionCluster aclusters[TrackerTraits::maxNumClustersPerModules];
0112 #ifdef EDM_ML_DEBUG
0113   auto totClustersFilled = 0;
0114 #endif
0115
0116   auto fillClusters = [&](uint32_t detId) {
0117     if (nclus < 0)
0118       return;  // this in reality should never happen
0119     edmNew::DetSetVector<SiPixelCluster>::FastFiller spc(*outputClusters, detId);
0120     auto layer = (DetId(detId).subdetId() == 1) ? ttopo.pxbLayer(detId) : 0;
0121     auto clusterThreshold = clusterThresholds_.getThresholdForLayerOnCondition(layer == 1);
0122     for (int32_t ic = 0; ic < nclus + 1; ++ic) {
0123       auto const& acluster = aclusters[ic];
0124       // in any case we cannot  go out of sync with gpu...
0125       if (!std::is_base_of<pixelTopology::Phase2, TrackerTraits>::value and acluster.charge < clusterThreshold)
0126         edm::LogWarning("SiPixelDigisClustersFromSoA") << "cluster below charge Threshold "
0127                                                        << "Layer/DetId/clusId " << layer << '/' << detId << '/' << ic
0128                                                        << " size/charge " << acluster.isize << '/' << acluster.charge;
0129       // sort by row (x)
0130       spc.emplace_back(acluster.isize, acluster.adc, acluster.x, acluster.y, acluster.xmin, acluster.ymin, ic);
0131       aclusters[ic].clear();
0132 #ifdef EDM_ML_DEBUG
0133       ++totClustersFilled;
0134       const auto& cluster{spc.back()};
0135       LogDebug("SiPixelDigisClustersFromSoA")
0136           << "putting in this cluster " << ic << " " << cluster.charge() << " " << cluster.pixelADC().size();
0137 #endif
0138       std::push_heap(spc.begin(), spc.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0139         return cl1.minPixelRow() < cl2.minPixelRow();
0140       });
0141     }
0142     nclus = -1;
0143     // sort by row (x)
0144     std::sort_heap(spc.begin(), spc.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0145       return cl1.minPixelRow() < cl2.minPixelRow();
0146     });
0147     if (spc.empty())
0148       spc.abort();
0149   };
0150
0151   for (uint32_t i = 0; i < nDigis; i++) {
0152     // check for uninitialized digis
0153     if (digis.rawIdArr(i) == 0)
0154       continue;
0155     // check for noisy/dead pixels (electrons set to 0)
0156     if (digis.adc(i) == 0)
0157       continue;
0158     if (digis.clus(i) > 9000)
0159       continue;  // not in cluster; TODO add an assert for the size
0160 #ifdef EDM_ML_DEBUG
0161     assert(digis.rawIdArr(i) > 109999);
0162 #endif
0163     if (detId != digis.rawIdArr(i)) {
0164       // new module
0165       fillClusters(detId);
0166 #ifdef EDM_ML_DEBUG
0167       assert(nclus == -1);
0168 #endif
0169       detId = digis.rawIdArr(i);
0170       if (storeDigis_) {
0171         detDigis = &collection->find_or_insert(detId);
0172         if ((*detDigis).empty())
0173           (*detDigis).data.reserve(64);  // avoid the first relocations
0174         else {
0175           edm::LogWarning("SiPixelDigisClustersFromSoA")
0176               << "Problem det present twice in input! " << (*detDigis).detId();
0177         }
0178       }
0179     }
0180     PixelDigi dig(digis.pdigi(i));
0181     if (storeDigis_)
0182       (*detDigis).data.emplace_back(dig);
0183       // fill clusters
0184 #ifdef EDM_ML_DEBUG
0185     assert(digis.clus(i) >= 0);
0186     assert(digis.clus(i) < static_cast<int32_t>(TrackerTraits::maxNumClustersPerModules));
0187 #endif
0188     nclus = std::max(digis.clus(i), nclus);
0189     auto row = dig.row();
0190     auto col = dig.column();
0191     SiPixelCluster::PixelPos pix(row, col);
0192     aclusters[digis.clus(i)].add(pix, digis.adc(i));
0193   }
0194
0195   // fill final clusters
0196   if (detId > 0)
0197     fillClusters(detId);
0198
0199 #ifdef EDM_ML_DEBUG
0200   LogDebug("SiPixelDigisClustersFromSoA") << "filled " << totClustersFilled << " clusters";
0201 #endif
0202
0203   if (produceDigis_)
0204     iEvent.put(digiPutToken_, std::move(collection));
0205   iEvent.put(clusterPutToken_, std::move(outputClusters));
0206 }
0207
0208 using SiPixelDigisClustersFromSoAPhase1 = SiPixelDigisClustersFromSoAT<pixelTopology::Phase1>;
0209 DEFINE_FWK_MODULE(SiPixelDigisClustersFromSoAPhase1);
0210 using SiPixelDigisClustersFromSoAPhase2 = SiPixelDigisClustersFromSoAT<pixelTopology::Phase2>;
0211 DEFINE_FWK_MODULE(SiPixelDigisClustersFromSoAPhase2);
0212 using SiPixelDigisClustersFromSoAHIonPhase1 = SiPixelDigisClustersFromSoAT<pixelTopology::HIonPhase1>;
0213 DEFINE_FWK_MODULE(SiPixelDigisClustersFromSoAHIonPhase1);