plugins/alpaka/HGCalSoARecHitsProducer.cc

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