HGCalRecProducers/plugins/HGCalLayerClustersFromSoAProducer.cc

0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002
0003 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0004
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0009 #include "FWCore/ParameterSet/interface/allowedValues.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/Utilities/interface/EDGetToken.h"
0012
0013 #include "DataFormats/Common/interface/ValueMap.h"
0014 #include "DataFormats/Math/interface/Point3D.h"
0015 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0016 #include "DataFormats/HGCalReco/interface/HGCalSoAClustersHostCollection.h"
0017 #include "DataFormats/HGCalReco/interface/HGCalSoARecHitsExtraHostCollection.h"
0018 #include "DataFormats/HGCalReco/interface/HGCalSoARecHitsHostCollection.h"
0019
0020 #include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"
0021
0022 #define DEBUG_CLUSTERS_ALPAKA 0
0023
0024 #if DEBUG_CLUSTERS_ALPAKA
0025 #include "RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h"
0026 #endif
0027
0028 class HGCalLayerClustersFromSoAProducer : public edm::stream::EDProducer<> {
0029 public:
0030   HGCalLayerClustersFromSoAProducer(edm::ParameterSet const& config)
0031       : getTokenSoAClusters_(consumes(config.getParameter<edm::InputTag>("src"))),
0032         getTokenSoACells_(consumes(config.getParameter<edm::InputTag>("hgcalRecHitsSoA"))),
0033         getTokenSoARecHitsExtra_(consumes(config.getParameter<edm::InputTag>("hgcalRecHitsLayerClustersSoA"))),
0034         detector_(config.getParameter<std::string>("detector")),
0035         hitsTime_(config.getParameter<unsigned int>("nHitsTime")),
0036         timeClname_(config.getParameter<std::string>("timeClname")) {
0037 #if DEBUG_CLUSTERS_ALPAKA
0038     moduleLabel_ = config.getParameter<std::string>("@module_label");
0039 #endif
0040     if (detector_ == "HFNose") {
0041       algoId_ = reco::CaloCluster::hfnose;
0042     } else if (detector_ == "EE") {
0043       algoId_ = reco::CaloCluster::hgcal_em;
0044     } else {  //for FH or BH
0045       algoId_ = reco::CaloCluster::hgcal_had;
0046     }
0047
0048     produces<std::vector<float>>("InitialLayerClustersMask");
0049     produces<std::vector<reco::BasicCluster>>();
0050     produces<edm::ValueMap<std::pair<float, float>>>(timeClname_);
0051   }
0052
0053   ~HGCalLayerClustersFromSoAProducer() override = default;
0054
0055   void produce(edm::Event& iEvent, edm::EventSetup const& iSetup) override {
0056     auto const& deviceData = iEvent.get(getTokenSoAClusters_);
0057
0058     auto const& deviceSoARecHitsExtra = iEvent.get(getTokenSoARecHitsExtra_);
0059     auto const soaRecHitsExtra_v = deviceSoARecHitsExtra.view();
0060
0061     auto const& deviceSoACells = iEvent.get(getTokenSoACells_);
0062     auto const soaCells_v = deviceSoACells.view();
0063
0064     auto const deviceView = deviceData.view();
0065
0066     std::unique_ptr<std::vector<reco::BasicCluster>> clusters(new std::vector<reco::BasicCluster>);
0067     clusters->reserve(deviceData->metadata().size());
0068
0069     // Create a vector of <clusters> locations, where each location holds a
0070     // vector of <nCells> floats. These vectors are used to compute the time for
0071     // each cluster.
0072     std::vector<std::vector<float>> times(deviceData->metadata().size());
0073     std::vector<std::vector<float>> timeErrors(deviceData->metadata().size());
0074
0075     for (int i = 0; i < deviceData->metadata().size(); ++i) {
0076       std::vector<std::pair<DetId, float>> thisCluster;
0077       thisCluster.reserve(deviceView.cells(i));
0078       clusters->emplace_back(deviceView.energy(i),
0079                              math::XYZPoint(deviceView.x(i), deviceView.y(i), deviceView.z(i)),
0080                              reco::CaloID::DET_HGCAL_ENDCAP,
0081                              std::move(thisCluster),
0082                              algoId_);
0083       clusters->back().setSeed(deviceView.seed(i));
0084       times[i].reserve(deviceView.cells(i));
0085       timeErrors[i].reserve(deviceView.cells(i));
0086     }
0087
0088     // Populate hits and fractions required to compute the cluster's time.
0089     // This procedure is complex and involves two SoAs: the original RecHits
0090     // SoA and the clustering algorithm's output SoA. Both SoAs have the same
0091     // cardinality, and crucially, the output SoA includes the cluster index.
0092     for (int32_t i = 0; i < soaRecHitsExtra_v.metadata().size(); ++i) {
0093       if (soaRecHitsExtra_v[i].clusterIndex() == -1) {
0094         continue;
0095       }
0096       assert(soaRecHitsExtra_v[i].clusterIndex() < (int)clusters->size());
0097       (*clusters)[soaRecHitsExtra_v[i].clusterIndex()].addHitAndFraction(soaCells_v[i].detid(), 1.f);
0098       if (soaCells_v[i].timeError() < 0.f) {
0099         continue;
0100       }
0101       times[soaRecHitsExtra_v[i].clusterIndex()].push_back(soaCells_v[i].time());
0102       timeErrors[soaRecHitsExtra_v[i].clusterIndex()].push_back(
0103           1.f / (soaCells_v[i].timeError() * soaCells_v[i].timeError()));
0104     }
0105
0106     // Finally, compute and assign the time to each cluster.
0107     std::vector<std::pair<float, float>> cluster_times;
0108     cluster_times.reserve(clusters->size());
0109     hgcalsimclustertime::ComputeClusterTime timeEstimator;
0110     for (unsigned i = 0; i < clusters->size(); ++i) {
0111       if (detector_ != "BH") {
0112         cluster_times.push_back(timeEstimator.fixSizeHighestDensity(times[i], timeErrors[i], hitsTime_));
0113       } else {
0114         cluster_times.push_back(std::pair<float, float>(-99.f, -1.f));
0115       }
0116     }
0117
0118 #if DEBUG_CLUSTERS_ALPAKA
0119     auto runNumber = iEvent.eventAuxiliary().run();
0120     auto lumiNumber = iEvent.eventAuxiliary().luminosityBlock();
0121     auto evtNumber = iEvent.eventAuxiliary().id().event();
0122
0123     hgcalUtils::DumpCellsSoA dumperCellsSoA;
0124     dumperCellsSoA.dumpInfos(deviceSoACells, moduleLabel_, runNumber, lumiNumber, evtNumber);
0125
0126     hgcalUtils::DumpClusters dumper;
0127     dumper.dumpInfos(*clusters, moduleLabel_, runNumber, lumiNumber, evtNumber, true);
0128
0129     hgcalUtils::DumpClustersSoA dumperSoA;
0130     dumperSoA.dumpInfos(deviceSoARecHitsExtra, moduleLabel_, runNumber, lumiNumber, evtNumber);
0131 #endif
0132
0133     auto clusterHandle = iEvent.put(std::move(clusters));
0134
0135     auto timeCl = std::make_unique<edm::ValueMap<std::pair<float, float>>>();
0136     edm::ValueMap<std::pair<float, float>>::Filler filler(*timeCl);
0137     filler.insert(clusterHandle, cluster_times.begin(), cluster_times.end());
0138     filler.fill();
0139     iEvent.put(std::move(timeCl), timeClname_);
0140
0141     // The layerClusterMask for the HGCAL detector is created at a later
0142     // stage, when the layer clusters from the different components of HGCAL
0143     // are merged together into a unique collection. For the case of HFNose,
0144     // since there is no further merging step needed, we create the
0145     // layerClustersMask directly here.
0146     if (detector_ == "HFNose") {
0147       std::unique_ptr<std::vector<float>> layerClustersMask(new std::vector<float>);
0148       layerClustersMask->resize(clusterHandle->size(), 1.0);
0149       iEvent.put(std::move(layerClustersMask), "InitialLayerClustersMask");
0150     }
0151   }
0152
0153   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0154     edm::ParameterSetDescription desc;
0155     desc.add<edm::InputTag>("src", edm::InputTag("hltHgcalSoALayerClustersProducer"));
0156     desc.add<edm::InputTag>("hgcalRecHitsLayerClustersSoA", edm::InputTag("hltHgcalSoARecHitsLayerClustersProducer"));
0157     desc.add<edm::InputTag>("hgcalRecHitsSoA", edm::InputTag("hltHgcalSoARecHitsProducer"));
0158     desc.add<unsigned int>("nHitsTime", 3);
0159     desc.add<std::string>("timeClname", "timeLayerCluster");
0160     desc.ifValue(edm::ParameterDescription<std::string>(
0161                      "detector", "EE", true, edm::Comment("the HGCAL component used to create clusters.")),
0162                  edm::allowedValues<std::string>("EE", "FH"));
0163     descriptions.addWithDefaultLabel(desc);
0164   }
0165
0166 private:
0167   edm::EDGetTokenT<HGCalSoAClustersHostCollection> const getTokenSoAClusters_;
0168   edm::EDGetTokenT<HGCalSoARecHitsHostCollection> const getTokenSoACells_;
0169   edm::EDGetTokenT<HGCalSoARecHitsExtraHostCollection> const getTokenSoARecHitsExtra_;
0170   std::string detector_;
0171   unsigned int hitsTime_;
0172   std::string timeClname_;
0173   reco::CaloCluster::AlgoId algoId_;
0174 #if DEBUG_CLUSTERS_ALPAKA
0175   std::string moduleLabel_;
0176 #endif
0177 };
0178 #include "FWCore/Framework/interface/MakerMacros.h"
0179 DEFINE_FWK_MODULE(HGCalLayerClustersFromSoAProducer);