plugins/alpaka/HGCalSoARecHitsProducer.cc

0001 #include "DataFormats/PortableTestObjects/interface/alpaka/TestDeviceCollection.h"
0002 #include "FWCore/Framework/interface/ConsumesCollector.h"
0003 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/EDProducer.h"
0008 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/EDPutToken.h"
0009 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h"
0010 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0011 #include "HeterogeneousCore/AlpakaTest/interface/AlpakaESTestRecords.h"
0012 #include "HeterogeneousCore/AlpakaTest/interface/alpaka/AlpakaESTestData.h"
0013 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0014 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0015 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0016
0017 #include "DataFormats/HGCalReco/interface/HGCalSoARecHitsHostCollection.h"
0018 #include "DataFormats/HGCalReco/interface/alpaka/HGCalSoARecHitsDeviceCollection.h"
0019
0020 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0021
0022   class HGCalSoARecHitsProducer : public stream::EDProducer<> {
0023   public:
0024     HGCalSoARecHitsProducer(edm::ParameterSet const& config)
0025         : detector_(config.getParameter<std::string>("detector")),
0026           initialized_(false),
0027           isNose_(detector_ == "HFNose"),
0028           maxNumberOfThickIndices_(config.getParameter<unsigned>("maxNumberOfThickIndices")),
0029           fcPerEle_(config.getParameter<double>("fcPerEle")),
0030           ecut_(config.getParameter<double>("ecut")),
0031           fcPerMip_(config.getParameter<std::vector<double>>("fcPerMip")),
0032           nonAgedNoises_(config.getParameter<std::vector<double>>("noises")),
0033           dEdXweights_(config.getParameter<std::vector<double>>("dEdXweights")),
0034           thicknessCorrection_(config.getParameter<std::vector<double>>("thicknessCorrection")),
0035           caloGeomToken_(consumesCollector().esConsumes<CaloGeometry, CaloGeometryRecord>()),
0036           hits_token_(consumes<HGCRecHitCollection>(config.getParameter<edm::InputTag>("recHits"))),
0037           deviceToken_{produces()} {}
0038
0039     ~HGCalSoARecHitsProducer() override = default;
0040
0041     void produce(device::Event& iEvent, device::EventSetup const& iSetup) override {
0042       edm::Handle<HGCRecHitCollection> hits_h;
0043
0044       edm::ESHandle<CaloGeometry> geom = iSetup.getHandle(caloGeomToken_);
0045       rhtools_.setGeometry(*geom);
0046       maxlayer_ = rhtools_.lastLayer(isNose_);
0047
0048       hits_h = iEvent.getHandle(hits_token_);
0049       auto const& hits = *(hits_h.product());
0050       computeThreshold();
0051
0052       // Count effective hits above threshold
0053       uint32_t index = 0;
0054       for (unsigned int i = 0; i < hits.size(); ++i) {
0055         const HGCRecHit& hgrh = hits[i];
0056         DetId detid = hgrh.detid();
0057         unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
0058
0059         // set sigmaNoise default value 1 to use kappa value directly in case of
0060         // sensor-independent thresholds
0061         int thickness_index = rhtools_.getSiThickIndex(detid);
0062         if (thickness_index == -1) {
0063           thickness_index = maxNumberOfThickIndices_;
0064         }
0065         double storedThreshold = thresholds_[layerOnSide][thickness_index];
0066         if (hgrh.energy() < storedThreshold)
0067           continue;  // this sets the ZS threshold at ecut times the sigma noise
0068         index++;
0069       }
0070
0071       // Allocate Host SoA will contain one entry for each RecHit above threshold
0072       HGCalSoARecHitsHostCollection cells(index, iEvent.queue());
0073       auto cellsView = cells.view();
0074
0075       // loop over all hits and create the Hexel structure, skip energies below ecut
0076       // for each layer and wafer calculate the thresholds (sigmaNoise and energy)
0077       // once
0078       index = 0;
0079       for (unsigned int i = 0; i < hits.size(); ++i) {
0080         const HGCRecHit& hgrh = hits[i];
0081         DetId detid = hgrh.detid();
0082         unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
0083
0084         // set sigmaNoise default value 1 to use kappa value directly in case of
0085         // sensor-independent thresholds
0086         float sigmaNoise = 1.f;
0087         int thickness_index = rhtools_.getSiThickIndex(detid);
0088         if (thickness_index == -1) {
0089           thickness_index = maxNumberOfThickIndices_;
0090         }
0091         double storedThreshold = thresholds_[layerOnSide][thickness_index];
0092         if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) {
0093           storedThreshold = thresholds_.at(layerOnSide).at(thickness_index + deltasi_index_regemfac_);
0094         }
0095         sigmaNoise = v_sigmaNoise_.at(layerOnSide).at(thickness_index);
0096
0097         if (hgrh.energy() < storedThreshold)
0098           continue;  // this sets the ZS threshold at ecut times the sigma noise
0099         // for the sensor
0100
0101         const GlobalPoint position(rhtools_.getPosition(detid));
0102         int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_;
0103         int layer = layerOnSide + offset;
0104         auto entryInSoA = cellsView[index];
0105         if (detector_ == "BH") {
0106           entryInSoA.dim1() = position.eta();
0107           entryInSoA.dim2() = position.phi();
0108         }  // else, isSilicon == true and eta phi values will not be used
0109         else {
0110           entryInSoA.dim1() = position.x();
0111           entryInSoA.dim2() = position.y();
0112         }
0113         entryInSoA.dim3() = position.z();
0114         entryInSoA.weight() = hgrh.energy();
0115         entryInSoA.sigmaNoise() = sigmaNoise;
0116         entryInSoA.layer() = layer;
0117         entryInSoA.recHitIndex() = i;
0118         entryInSoA.detid() = detid.rawId();
0119         entryInSoA.time() = hgrh.time();
0120         entryInSoA.timeError() = hgrh.timeError();
0121         index++;
0122       }
0123 #if 0
0124         std::cout << "Size: " << cells->metadata().size() << " count cells: " << index
0125           << " i.e. " << cells->metadata().size() << std::endl;
0126 #endif
0127
0128       if constexpr (!std::is_same_v<Device, alpaka_common::DevHost>) {
0129         // Trigger copy async to GPU
0130         //std::cout << "GPU" << std::endl;
0131         HGCalSoARecHitsDeviceCollection deviceProduct{cells->metadata().size(), iEvent.queue()};
0132         alpaka::memcpy(iEvent.queue(), deviceProduct.buffer(), cells.const_buffer());
0133         iEvent.emplace(deviceToken_, std::move(deviceProduct));
0134       } else {
0135         //std::cout << "CPU" << std::endl;
0136         iEvent.emplace(deviceToken_, std::move(cells));
0137       }
0138     }
0139
0140     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0141       edm::ParameterSetDescription desc;
0142       desc.add<std::string>("detector", "EE")->setComment("options EE, FH, BH,  HFNose; other value defaults to EE");
0143       desc.add<edm::InputTag>("recHits", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0144       desc.add<unsigned int>("maxNumberOfThickIndices", 6);
0145       desc.add<double>("fcPerEle", 0.00016020506);
0146       desc.add<std::vector<double>>("fcPerMip");
0147       desc.add<std::vector<double>>("thicknessCorrection");
0148       desc.add<std::vector<double>>("noises");
0149       desc.add<std::vector<double>>("dEdXweights");
0150       desc.add<double>("ecut", 3.);
0151       descriptions.addWithDefaultLabel(desc);
0152     }
0153
0154   private:
0155     std::string detector_;
0156     bool initialized_;
0157     bool isNose_;
0158     unsigned maxNumberOfThickIndices_;
0159     unsigned int maxlayer_;
0160     int deltasi_index_regemfac_;
0161     double sciThicknessCorrection_;
0162     double fcPerEle_;
0163     double ecut_;
0164     std::vector<double> fcPerMip_;
0165     std::vector<double> nonAgedNoises_;
0166     std::vector<double> dEdXweights_;
0167     std::vector<double> thicknessCorrection_;
0168     std::vector<std::vector<double>> thresholds_;
0169     std::vector<std::vector<double>> v_sigmaNoise_;
0170
0171     hgcal::RecHitTools rhtools_;
0172     edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0173     edm::EDGetTokenT<HGCRecHitCollection> hits_token_;
0174     device::EDPutToken<HGCalSoARecHitsDeviceCollection> const deviceToken_;
0175
0176     void computeThreshold() {
0177       // To support the TDR geometry and also the post-TDR one (v9 onwards), we
0178       // need to change the logic of the vectors containing signal to noise and
0179       // thresholds. The first 3 indices will keep on addressing the different
0180       // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address
0181       // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the
0182       // seventh) will address the Scintillators. This change will support both
0183       // geometries at the same time.
0184
0185       if (initialized_)
0186         return;  // only need to calculate thresholds once
0187
0188       initialized_ = true;
0189
0190       std::vector<double> dummy;
0191
0192       dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0);  // +1 to accomodate for the Scintillators
0193       thresholds_.resize(maxlayer_, dummy);
0194       v_sigmaNoise_.resize(maxlayer_, dummy);
0195
0196       for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
0197         for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
0198           float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
0199                              (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
0200           thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
0201           v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
0202 #if 0
0203             std::cout << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick]
0204               << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick]
0205               << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick]
0206               << " sigmaNoise: " << sigmaNoise << "\n";
0207 #endif
0208         }
0209       }
0210     }
0211   };
0212
0213 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE
0214
0215 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h"
0216 DEFINE_FWK_ALPAKA_MODULE(HGCalSoARecHitsProducer);