SiPixelRecHits/plugins/SiPixelRecHitSoAFromLegacy.cc

0001 #include <cuda_runtime.h>
0002
0003 #include "CUDADataFormats/BeamSpot/interface/BeamSpotCUDA.h"
0004 #include "CUDADataFormats/SiPixelCluster/interface/SiPixelClustersCUDA.h"
0005 #include "CUDADataFormats/SiPixelDigi/interface/SiPixelDigisCUDA.h"
0006 #include "CUDADataFormats/TrackingRecHit/interface/TrackingRecHitSoAHost.h"
0007 #include "CUDADataFormats/Common/interface/PortableHostCollection.h"
0008 #include "CUDADataFormats/Common/interface/HostProduct.h"
0009 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0010 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/Framework/interface/global/EDProducer.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0023 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0024 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0025 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0026 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h"
0027 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFast.h"
0028
0029 #include "gpuPixelRecHits.h"
0030
0031 template <typename TrackerTraits>
0032 class SiPixelRecHitSoAFromLegacyT : public edm::global::EDProducer<> {
0033 public:
0034   explicit SiPixelRecHitSoAFromLegacyT(const edm::ParameterSet& iConfig);
0035   ~SiPixelRecHitSoAFromLegacyT() override = default;
0036
0037   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0038
0039   using HitModuleStart = std::array<uint32_t, TrackerTraits::numberOfModules + 1>;
0040   using HMSstorage = HostProduct<uint32_t[]>;
0041   using HitsOnHost = TrackingRecHitSoAHost<TrackerTraits>;
0042
0043 private:
0044   void produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0045
0046   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0047   const edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> cpeToken_;
0048   const edm::EDGetTokenT<reco::BeamSpot> bsGetToken_;
0049   const edm::EDGetTokenT<SiPixelClusterCollectionNew> clusterToken_;  // Legacy Clusters
0050   const edm::EDPutTokenT<HitsOnHost> tokenHit_;
0051   const edm::EDPutTokenT<HMSstorage> tokenModuleStart_;
0052   const bool convert2Legacy_;
0053 };
0054
0055 template <typename TrackerTraits>
0056 SiPixelRecHitSoAFromLegacyT<TrackerTraits>::SiPixelRecHitSoAFromLegacyT(const edm::ParameterSet& iConfig)
0057     : geomToken_(esConsumes()),
0058       cpeToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("CPE")))),
0059       bsGetToken_{consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))},
0060       clusterToken_{consumes<SiPixelClusterCollectionNew>(iConfig.getParameter<edm::InputTag>("src"))},
0061       tokenHit_{produces<HitsOnHost>()},
0062       tokenModuleStart_{produces<HMSstorage>()},
0063       convert2Legacy_(iConfig.getParameter<bool>("convertToLegacy")) {
0064   if (convert2Legacy_)
0065     produces<SiPixelRecHitCollectionNew>();
0066 }
0067
0068 template <typename TrackerTraits>
0069 void SiPixelRecHitSoAFromLegacyT<TrackerTraits>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0070   edm::ParameterSetDescription desc;
0071
0072   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0073   desc.add<edm::InputTag>("src", edm::InputTag("siPixelClustersPreSplitting"));
0074   std::string cpeName = "PixelCPEFast";
0075   cpeName += TrackerTraits::nameModifier;
0076   desc.add<std::string>("CPE", cpeName);
0077   desc.add<bool>("convertToLegacy", false);
0078
0079   descriptions.addWithDefaultLabel(desc);
0080 }
0081
0082 template <typename TrackerTraits>
0083 void SiPixelRecHitSoAFromLegacyT<TrackerTraits>::produce(edm::StreamID streamID,
0084                                                          edm::Event& iEvent,
0085                                                          const edm::EventSetup& es) const {
0086   const TrackerGeometry* geom_ = &es.getData(geomToken_);
0087   PixelCPEFast<TrackerTraits> const* fcpe = dynamic_cast<const PixelCPEFast<TrackerTraits>*>(&es.getData(cpeToken_));
0088   if (not fcpe) {
0089     throw cms::Exception("Configuration") << "SiPixelRecHitSoAFromLegacy can only use a CPE of type PixelCPEFast";
0090   }
0091   auto const& cpeView = fcpe->getCPUProduct();
0092
0093   const reco::BeamSpot& bs = iEvent.get(bsGetToken_);
0094
0095   BeamSpotPOD bsHost;
0096   bsHost.x = bs.x0();
0097   bsHost.y = bs.y0();
0098   bsHost.z = bs.z0();
0099
0100   edm::Handle<SiPixelClusterCollectionNew> hclusters;
0101   iEvent.getByToken(clusterToken_, hclusters);
0102   auto const& input = *hclusters;
0103
0104   constexpr int nModules = TrackerTraits::numberOfModules;
0105   constexpr int startBPIX2 = pixelTopology::layerStart<TrackerTraits>(1);
0106
0107   // allocate a buffer for the indices of the clusters
0108   auto hmsp = std::make_unique<uint32_t[]>(nModules + 1);
0109   auto hitsModuleStart = hmsp.get();
0110   // wrap the buffer in a HostProduct
0111   auto hms = std::make_unique<HMSstorage>(std::move(hmsp));
0112   // move the HostProduct to the Event, without reallocating the buffer or affecting hitsModuleStart
0113   iEvent.put(tokenModuleStart_, std::move(hms));
0114
0115   // legacy output
0116   auto legacyOutput = std::make_unique<SiPixelRecHitCollectionNew>();
0117
0118   std::vector<edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster>> clusterRef;
0119
0120   constexpr uint32_t maxHitsInModule = TrackerTraits::maxNumClustersPerModules;
0121
0122   cms::cuda::PortableHostCollection<SiPixelClustersCUDALayout<>> clusters_h(nModules + 1);
0123
0124   memset(clusters_h.view().clusInModule(), 0, (nModules + 1) * sizeof(uint32_t));  // needed??
0125   memset(clusters_h.view().moduleStart(), 0, (nModules + 1) * sizeof(uint32_t));
0126   memset(clusters_h.view().moduleId(), 0, (nModules + 1) * sizeof(uint32_t));
0127   memset(clusters_h.view().clusModuleStart(), 0, (nModules + 1) * sizeof(uint32_t));
0128
0129   assert(0 == clusters_h.view()[nModules].clusInModule());
0130   clusters_h.view()[1].moduleStart() = 0;
0131
0132   // fill cluster arrays
0133   int numberOfClusters = 0;
0134   for (auto const& dsv : input) {
0135     unsigned int detid = dsv.detId();
0136     DetId detIdObject(detid);
0137     const GeomDetUnit* genericDet = geom_->idToDetUnit(detIdObject);
0138     auto gind = genericDet->index();
0139     assert(gind < nModules);
0140     auto const nclus = dsv.size();
0141     clusters_h.view()[gind].clusInModule() = nclus;
0142     numberOfClusters += nclus;
0143   }
0144   clusters_h.view()[0].clusModuleStart() = 0;
0145
0146   for (int i = 1; i < nModules + 1; ++i) {
0147     clusters_h.view()[i].clusModuleStart() =
0148         clusters_h.view()[i - 1].clusModuleStart() + clusters_h.view()[i - 1].clusInModule();
0149   }
0150
0151   assert((uint32_t)numberOfClusters == clusters_h.view()[nModules].clusModuleStart());
0152   // output SoA
0153   // element 96 is the start of BPIX2 (i.e. the number of clusters in BPIX1)
0154   HitsOnHost output(
0155       numberOfClusters, clusters_h.view()[startBPIX2].clusModuleStart(), &cpeView, clusters_h.view().clusModuleStart());
0156
0157   if (0 == numberOfClusters) {
0158     iEvent.emplace(tokenHit_, std::move(output));
0159     if (convert2Legacy_)
0160       iEvent.put(std::move(legacyOutput));
0161     return;
0162   }
0163
0164   if (convert2Legacy_)
0165     legacyOutput->reserve(nModules, numberOfClusters);
0166
0167   int numberOfDetUnits = 0;
0168   int numberOfHits = 0;
0169   for (auto const& dsv : input) {
0170     numberOfDetUnits++;
0171     unsigned int detid = dsv.detId();
0172     DetId detIdObject(detid);
0173     const GeomDetUnit* genericDet = geom_->idToDetUnit(detIdObject);
0174     auto const gind = genericDet->index();
0175     assert(gind < nModules);
0176     const PixelGeomDetUnit* pixDet = dynamic_cast<const PixelGeomDetUnit*>(genericDet);
0177     assert(pixDet);
0178     auto const nclus = dsv.size();
0179
0180     assert(clusters_h.view()[gind].clusInModule() == nclus);
0181     if (0 == nclus)
0182       continue;  // is this really possible?
0183
0184     auto const fc = clusters_h.view()[gind].clusModuleStart();
0185     auto const lc = clusters_h.view()[gind + 1].clusModuleStart();
0186     assert(lc > fc);
0187     LogDebug("SiPixelRecHitSoAFromLegacy") << "in det " << gind << ": conv " << nclus << " hits from " << dsv.size()
0188                                            << " legacy clusters" << ' ' << fc << ',' << lc;
0189     assert((lc - fc) == nclus);
0190     if (nclus > maxHitsInModule)
0191       printf(
0192           "WARNING: too many clusters %d in Module %d. Only first %d Hits converted\n", nclus, gind, maxHitsInModule);
0193
0194     // count digis
0195     uint32_t ndigi = 0;
0196     for (auto const& clust : dsv) {
0197       assert(clust.size() > 0);
0198       ndigi += clust.size();
0199     }
0200
0201     cms::cuda::PortableHostCollection<SiPixelDigisSoA> digis_h(ndigi);
0202
0203     clusterRef.clear();
0204     clusters_h.view()[0].moduleId() = gind;
0205
0206     uint32_t ic = 0;
0207     ndigi = 0;
0208     //filling digis
0209     for (auto const& clust : dsv) {
0210       assert(clust.size() > 0);
0211       for (int i = 0, nd = clust.size(); i < nd; ++i) {
0212         auto px = clust.pixel(i);
0213         digis_h.view()[ndigi].xx() = px.x;
0214         digis_h.view()[ndigi].yy() = px.y;
0215         digis_h.view()[ndigi].adc() = px.adc;
0216         digis_h.view()[ndigi].moduleId() = gind;
0217         digis_h.view()[ndigi].clus() = ic;
0218         ++ndigi;
0219       }
0220
0221       if (convert2Legacy_)
0222         clusterRef.emplace_back(edmNew::makeRefTo(hclusters, &clust));
0223       ic++;
0224     }
0225     assert(nclus == ic);
0226
0227     numberOfHits += nclus;
0228     // filled creates view
0229     assert(digis_h.view()[0].adc() != 0);
0230     // we run on blockId.x==0
0231
0232     gpuPixelRecHits::getHits(&cpeView, &bsHost, digis_h.view(), ndigi, clusters_h.view(), output.view());
0233     for (auto h = fc; h < lc; ++h)
0234       if (h - fc < maxHitsInModule)
0235         assert(gind == output.view()[h].detectorIndex());
0236       else
0237         assert(gpuClustering::invalidModuleId == output.view()[h].detectorIndex());
0238
0239     if (convert2Legacy_) {
0240       SiPixelRecHitCollectionNew::FastFiller recHitsOnDetUnit(*legacyOutput, detid);
0241       for (auto h = fc; h < lc; ++h) {
0242         auto ih = h - fc;
0243
0244         if (ih >= maxHitsInModule)
0245           break;
0246
0247         assert(ih < clusterRef.size());
0248         LocalPoint lp(output.view()[h].xLocal(), output.view()[h].yLocal());
0249         LocalError le(output.view()[h].xerrLocal(), 0, output.view()[h].yerrLocal());
0250
0251         SiPixelRecHitQuality::QualWordType rqw = 0;
0252         SiPixelRecHit hit(lp, le, rqw, *genericDet, clusterRef[ih]);
0253         recHitsOnDetUnit.push_back(hit);
0254       }
0255     }
0256   }
0257
0258   assert(numberOfHits == numberOfClusters);
0259
0260   // fill data structure to support CA
0261   constexpr auto nLayers = TrackerTraits::numberOfLayers;
0262   for (auto i = 0U; i < nLayers + 1; ++i) {
0263     output.view().hitsLayerStart()[i] = clusters_h.view()[cpeView.layerGeometry().layerStart[i]].clusModuleStart();
0264     LogDebug("SiPixelRecHitSoAFromLegacy")
0265         << "Layer n." << i << " - starting at module: " << cpeView.layerGeometry().layerStart[i]
0266         << " - starts ad cluster: " << output->hitsLayerStart()[i] << "\n";
0267   }
0268
0269   cms::cuda::fillManyFromVector(&(output.view().phiBinner()),
0270                                 nLayers,
0271                                 output.view().iphi(),
0272                                 output.view().hitsLayerStart().data(),
0273                                 output.view().nHits(),
0274                                 256,
0275                                 output.view().phiBinnerStorage());
0276
0277   LogDebug("SiPixelRecHitSoAFromLegacy") << "created HitSoa for " << numberOfClusters << " clusters in "
0278                                          << numberOfDetUnits << " Dets"
0279                                          << "\n";
0280
0281   // copy pointer to data (SoA view) to allocated buffer
0282   memcpy(hitsModuleStart, clusters_h.view().clusModuleStart(), nModules * sizeof(uint32_t));
0283
0284   iEvent.emplace(tokenHit_, std::move(output));
0285   if (convert2Legacy_)
0286     iEvent.put(std::move(legacyOutput));
0287 }
0288
0289 using SiPixelRecHitSoAFromLegacyPhase1 = SiPixelRecHitSoAFromLegacyT<pixelTopology::Phase1>;
0290 DEFINE_FWK_MODULE(SiPixelRecHitSoAFromLegacyPhase1);
0291
0292 using SiPixelRecHitSoAFromLegacyPhase2 = SiPixelRecHitSoAFromLegacyT<pixelTopology::Phase2>;
0293 DEFINE_FWK_MODULE(SiPixelRecHitSoAFromLegacyPhase2);
0294
0295 using SiPixelRecHitSoAFromLegacyHIonPhase1 = SiPixelRecHitSoAFromLegacyT<pixelTopology::HIonPhase1>;
0296 DEFINE_FWK_MODULE(SiPixelRecHitSoAFromLegacyHIonPhase1);