diff --git a/DataFormats/ParticleFlowReco/interface/PFRecHitSoA.h b/DataFormats/ParticleFlowReco/interface/PFRecHitSoA.h index 19dd4c54fca8d..cdb12377e4019 100644 --- a/DataFormats/ParticleFlowReco/interface/PFRecHitSoA.h +++ b/DataFormats/ParticleFlowReco/interface/PFRecHitSoA.h @@ -14,6 +14,7 @@ namespace reco { using PFRecHitsNeighbours = Eigen::Matrix; GENERATE_SOA_LAYOUT(PFRecHitSoALayout, SOA_COLUMN(uint32_t, detId), + SOA_COLUMN(uint32_t, denseId), SOA_COLUMN(float, energy), SOA_COLUMN(float, time), SOA_COLUMN(int, depth), diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFClusterParamsSoA.h b/RecoParticleFlow/PFClusterProducer/interface/PFClusterParamsSoA.h index dded4580d0e1a..9003ce1c527b7 100644 --- a/RecoParticleFlow/PFClusterProducer/interface/PFClusterParamsSoA.h +++ b/RecoParticleFlow/PFClusterProducer/interface/PFClusterParamsSoA.h @@ -1,5 +1,5 @@ -#ifndef RecoParticleFlow_PFClusterProducer_interface_PFRecHitHBHEParamsSoA_h -#define RecoParticleFlow_PFClusterProducer_interface_PFRecHitHBHEParamsSoA_h +#ifndef RecoParticleFlow_PFClusterProducer_interface_PFClusterParamsSoA_h +#define RecoParticleFlow_PFClusterProducer_interface_PFClusterParamsSoA_h #include "DataFormats/SoATemplate/interface/SoACommon.h" #include "DataFormats/SoATemplate/interface/SoALayout.h" @@ -9,12 +9,12 @@ namespace reco { GENERATE_SOA_LAYOUT(PFClusterParamsSoALayout, SOA_SCALAR(int32_t, nNeigh), - SOA_SCALAR(float, seedPt2ThresholdEB), - SOA_SCALAR(float, seedPt2ThresholdEE), - SOA_COLUMN(float, seedEThresholdEB_vec), - SOA_COLUMN(float, seedEThresholdEE_vec), - SOA_COLUMN(float, topoEThresholdEB_vec), - SOA_COLUMN(float, topoEThresholdEE_vec), + SOA_SCALAR(float, seedPt2ThresholdHB), + SOA_SCALAR(float, seedPt2ThresholdHE), + SOA_COLUMN(float, seedEThresholdHB_vec), + SOA_COLUMN(float, seedEThresholdHE_vec), + SOA_COLUMN(float, topoEThresholdHB_vec), + SOA_COLUMN(float, topoEThresholdHE_vec), SOA_SCALAR(float, showerSigma2), SOA_SCALAR(float, minFracToKeep), SOA_SCALAR(float, minFracTot), @@ -23,8 +23,8 @@ namespace reco { SOA_SCALAR(float, stoppingTolerance), SOA_SCALAR(float, minFracInCalc), SOA_SCALAR(float, minAllowedNormalization), - SOA_COLUMN(float, recHitEnergyNormInvEB_vec), - SOA_COLUMN(float, recHitEnergyNormInvEE_vec), + SOA_COLUMN(float, recHitEnergyNormInvHB_vec), + SOA_COLUMN(float, recHitEnergyNormInvHE_vec), SOA_SCALAR(float, barrelTimeResConsts_corrTermLowE), SOA_SCALAR(float, barrelTimeResConsts_threshLowE), SOA_SCALAR(float, barrelTimeResConsts_noiseTerm), diff --git a/RecoParticleFlow/PFClusterProducer/plugins/LegacyPFClusterProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/LegacyPFClusterProducer.cc index 36204501d432d..4027bb340a168 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/LegacyPFClusterProducer.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/LegacyPFClusterProducer.cc @@ -22,6 +22,9 @@ #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Utilities/interface/EDPutToken.h" +#include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h" +#include "CondTools/Hcal/interface/HcalPFCutsHandler.h" + #include "DataFormats/Common/interface/Handle.h" #include "DataFormats/DetId/interface/DetId.h" #include "DataFormats/ParticleFlowReco/interface/PFRecHitHostCollection.h" @@ -39,7 +42,9 @@ class LegacyPFClusterProducer : public edm::stream::EDProducer<> { InputPFRecHitSoA_Token_{consumes(config.getParameter("PFRecHitsLabelIn"))}, pfClusParamsToken_(esConsumes(config.getParameter("pfClusterParams"))), legacyPfClustersToken_(produces()), - recHitsLabel_(consumes(config.getParameter("recHitsSource"))) { + recHitsLabel_(consumes(config.getParameter("recHitsSource"))), + hcalCutsToken_(esConsumes(edm::ESInputTag("", "withTopo"))), + cutsFromDB_(config.getParameter("usePFThresholdsFromDB")) { edm::ConsumesCollector cc = consumesCollector(); //setup pf cluster builder if requested @@ -65,6 +70,7 @@ class LegacyPFClusterProducer : public edm::stream::EDProducer<> { desc.add("PFRecHitsLabelIn"); desc.add("pfClusterParams"); desc.add("recHitsSource"); + desc.add("usePFThresholdsFromDB", true); { edm::ParameterSetDescription pfClusterBuilder; pfClusterBuilder.add("maxIterations", 5); @@ -180,6 +186,8 @@ class LegacyPFClusterProducer : public edm::stream::EDProducer<> { const edm::ESGetToken pfClusParamsToken_; const edm::EDPutTokenT legacyPfClustersToken_; const edm::EDGetTokenT recHitsLabel_; + const edm::ESGetToken hcalCutsToken_; + const bool cutsFromDB_; // the actual algorithm std::unique_ptr positionCalc_; std::unique_ptr allCellsPositionCalc_; @@ -188,6 +196,8 @@ class LegacyPFClusterProducer : public edm::stream::EDProducer<> { void LegacyPFClusterProducer::produce(edm::Event& event, const edm::EventSetup& setup) { const reco::PFRecHitHostCollection& pfRecHits = event.get(InputPFRecHitSoA_Token_); + HcalPFCuts const* paramPF = cutsFromDB_ ? &setup.getData(hcalCutsToken_) : nullptr; + auto const& pfClusterSoA = event.get(pfClusterSoAToken_).const_view(); auto const& pfRecHitFractionSoA = event.get(pfRecHitFractionSoAToken_).const_view(); @@ -221,11 +231,9 @@ void LegacyPFClusterProducer::produce(edm::Event& event, const edm::EventSetup& // Now PFRecHitFraction of this PFCluster is set. Now compute calculateAndSetPosition (energy, position etc) if (nTopoSeeds[pfClusterSoA[i].topoId()] == 1 && allCellsPositionCalc_) { - allCellsPositionCalc_->calculateAndSetPosition( - temp, nullptr); // temporarily use nullptr until we can properly set GT thresholds + allCellsPositionCalc_->calculateAndSetPosition(temp, paramPF); } else { - positionCalc_->calculateAndSetPosition( - temp, nullptr); // temporarily use nullptr until we can properly set GT thresholds + positionCalc_->calculateAndSetPosition(temp, paramPF); } out.emplace_back(std::move(temp)); } diff --git a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterParamsESProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterParamsESProducer.cc index 8a5da486752ce..36d4b025c3a48 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterParamsESProducer.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterParamsESProducer.cc @@ -35,17 +35,17 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { if (thresholds.size() != kMaxDepth_barrel) throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_barrel << ") for \"\" vector of det = \"" << det << "\""; - view.seedPt2ThresholdEB() = seedPt2Threshold; + view.seedPt2ThresholdHB() = seedPt2Threshold; for (size_t idx = 0; idx < thresholds.size(); ++idx) { - view.seedEThresholdEB_vec()[idx] = thresholds[idx]; + view.seedEThresholdHB_vec()[idx] = thresholds[idx]; } } else if (det == "HCAL_ENDCAP") { if (thresholds.size() != kMaxDepth_endcap) throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_endcap << ") for \"\" vector of det = \"" << det << "\""; - view.seedPt2ThresholdEE() = seedPt2Threshold; + view.seedPt2ThresholdHE() = seedPt2Threshold; for (size_t idx = 0; idx < thresholds.size(); ++idx) { - view.seedEThresholdEE_vec()[idx] = thresholds[idx]; + view.seedEThresholdHE_vec()[idx] = thresholds[idx]; } } else { throw cms::Exception("Configuration") << "Unknown detector when parsing seedFinder: " << det; @@ -63,14 +63,14 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_barrel << ") for \"\" vector of det = \"" << det << "\""; for (size_t idx = 0; idx < thresholds.size(); ++idx) { - view.topoEThresholdEB_vec()[idx] = thresholds[idx]; + view.topoEThresholdHB_vec()[idx] = thresholds[idx]; } } else if (det == "HCAL_ENDCAP") { if (thresholds.size() != kMaxDepth_endcap) throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_endcap << ") for \"\" vector of det = \"" << det << "\""; for (size_t idx = 0; idx < thresholds.size(); ++idx) { - view.topoEThresholdEE_vec()[idx] = thresholds[idx]; + view.topoEThresholdHE_vec()[idx] = thresholds[idx]; } } else { throw cms::Exception("Configuration") << "Unknown detector when parsing initClusteringStep: " << det; @@ -99,7 +99,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { << "Invalid size (" << recHitNorms.size() << " != " << kMaxDepth_barrel << ") for \"\" vector of det = \"" << det << "\""; for (size_t idx = 0; idx < recHitNorms.size(); ++idx) { - view.recHitEnergyNormInvEB_vec()[idx] = 1. / recHitNorms[idx]; + view.recHitEnergyNormInvHB_vec()[idx] = 1. / recHitNorms[idx]; } } else if (det == "HCAL_ENDCAP") { if (recHitNorms.size() != kMaxDepth_endcap) @@ -107,7 +107,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { << "Invalid size (" << recHitNorms.size() << " != " << kMaxDepth_endcap << ") for \"\" vector of det = \"" << det << "\""; for (size_t idx = 0; idx < recHitNorms.size(); ++idx) { - view.recHitEnergyNormInvEE_vec()[idx] = 1. / recHitNorms[idx]; + view.recHitEnergyNormInvHE_vec()[idx] = 1. / recHitNorms[idx]; } } else { throw cms::Exception("Configuration") << "Unknown detector when parsing recHitEnergyNorms: " << det; diff --git a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducer.cc index e4291f8e705c8..bde6db46b08d9 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducer.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducer.cc @@ -18,6 +18,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { public: PFClusterSoAProducer(edm::ParameterSet const& config) : pfClusParamsToken(esConsumes(config.getParameter("pfClusterParams"))), + topologyToken_(esConsumes(config.getParameter("topology"))), inputPFRecHitSoA_Token_{consumes(config.getParameter("pfRecHits"))}, outputPFClusterSoA_Token_{produces()}, outputPFRHFractionSoA_Token_{produces()}, @@ -26,6 +27,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { void produce(device::Event& event, device::EventSetup const& setup) override { const reco::PFClusterParamsDeviceCollection& params = setup.getData(pfClusParamsToken); + const reco::PFRecHitHCALTopologyDeviceCollection& topology = setup.getData(topologyToken_); const reco::PFRecHitHostCollection& pfRecHits = event.get(inputPFRecHitSoA_Token_); const int nRH = pfRecHits->size(); @@ -36,7 +38,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { PFClusterProducerKernel kernel(event.queue(), pfRecHits); kernel.execute( - event.queue(), params, pfClusteringVars, pfClusteringEdgeVars, pfRecHits, pfClusters, pfrhFractions); + event.queue(), params, topology, pfClusteringVars, pfClusteringEdgeVars, pfRecHits, pfClusters, pfrhFractions); if (synchronise_) alpaka::wait(event.queue()); @@ -49,6 +51,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { edm::ParameterSetDescription desc; desc.add("pfRecHits"); desc.add("pfClusterParams"); + desc.add("topology"); desc.add("synchronise"); desc.add("pfRecHitFractionAllocation", 120); descriptions.addWithDefaultLabel(desc); @@ -56,6 +59,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { private: const device::ESGetToken pfClusParamsToken; + const device::ESGetToken topologyToken_; const edm::EDGetTokenT inputPFRecHitSoA_Token_; const device::EDPutToken outputPFClusterSoA_Token_; const device::EDPutToken outputPFRHFractionSoA_Token_; diff --git a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.dev.cc b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.dev.cc index 6b0f4dc27d88b..ea7816ce0cb87 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.dev.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.dev.cc @@ -88,14 +88,16 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { // Processing single seed clusters // Device function designed to be called by all threads of a given block template >> - ALPAKA_FN_ACC static void hcalFastCluster_singleSeed(const TAcc& acc, - reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, - int topoId, // from selection - int nRHTopo, // from selection - reco::PFRecHitDeviceCollection::ConstView pfRecHits, - reco::PFClusteringVarsDeviceCollection::View pfClusteringVars, - reco::PFClusterDeviceCollection::View clusterView, - reco::PFRecHitFractionDeviceCollection::View fracView) { + ALPAKA_FN_ACC static void hcalFastCluster_singleSeed( + const TAcc& acc, + reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, + const reco::PFRecHitHCALTopologyDeviceCollection::ConstView topology, + int topoId, // from selection + int nRHTopo, // from selection + reco::PFRecHitDeviceCollection::ConstView pfRecHits, + reco::PFClusteringVarsDeviceCollection::View pfClusteringVars, + reco::PFClusterDeviceCollection::View clusterView, + reco::PFRecHitFractionDeviceCollection::View fracView) { int tid = alpaka::getIdx(acc)[0u]; // thread index is rechit number // Declaration of shared variables int& i = alpaka::declareSharedVar(acc); @@ -119,13 +121,17 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { clusterEnergy = seedEnergy; tol = pfClusParams.stoppingTolerance(); // stopping tolerance * tolerance scaling - if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) - rhENormInv = pfClusParams.recHitEnergyNormInvEB_vec()[pfRecHits[i].depth() - 1]; - else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) - rhENormInv = pfClusParams.recHitEnergyNormInvEE_vec()[pfRecHits[i].depth() - 1]; - else { - rhENormInv = 0.; - printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + if (topology.cutsFromDB()) { + rhENormInv = (1.f / topology[pfRecHits[i].denseId()].noiseThreshold()); + } else { + if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) + rhENormInv = pfClusParams.recHitEnergyNormInvHB_vec()[pfRecHits[i].depth() - 1]; + else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) + rhENormInv = pfClusParams.recHitEnergyNormInvHE_vec()[pfRecHits[i].depth() - 1]; + else { + rhENormInv = 0.; + printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + } } iter = 0; @@ -161,7 +167,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (clusterPos.z - rhPos.z) * (clusterPos.z - rhPos.z); d2 = dist2 / pfClusParams.showerSigma2(); - fraction = clusterEnergy * rhENormInv * expf(-0.5 * d2); + fraction = clusterEnergy * rhENormInv * expf(-0.5f * d2); // For single seed clusters, rechit fraction is either 1 (100%) or -1 (not included) if (fraction > pfClusParams.minFracTot() && d2 < cutoffDistance) @@ -253,6 +259,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ALPAKA_FN_ACC static void hcalFastCluster_multiSeedParallel( const TAcc& acc, reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, + const reco::PFRecHitHCALTopologyDeviceCollection::ConstView topology, int topoId, // from selection int nSeeds, // from selection int nRHTopo, // from selection @@ -289,12 +296,19 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { notDone = true; int i = pfClusteringVars[topoSeedBegin].topoSeedList(); - if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) - rhENormInv = pfClusParams.recHitEnergyNormInvEB_vec()[pfRecHits[i].depth() - 1]; - else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) - rhENormInv = pfClusParams.recHitEnergyNormInvEE_vec()[pfRecHits[i].depth() - 1]; - else - printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + + if (topology.cutsFromDB()) { + rhENormInv = (1.f / topology[pfRecHits[i].denseId()].noiseThreshold()); + } else { + if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) + rhENormInv = pfClusParams.recHitEnergyNormInvHB_vec()[pfRecHits[i].depth() - 1]; + else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) + rhENormInv = pfClusParams.recHitEnergyNormInvHE_vec()[pfRecHits[i].depth() - 1]; + else { + rhENormInv = 0.; + printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + } + } } alpaka::syncBlockThreads(acc); // all threads call sync @@ -391,7 +405,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (clusterPos[s].z - rhThreadPos.z) * (clusterPos[s].z - rhThreadPos.z); float d2 = dist2 / pfClusParams.showerSigma2(); - float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5 * d2); + float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5f * d2); rhFracSum[tid] += fraction; } @@ -406,7 +420,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (clusterPos[s].z - rhThreadPos.z) * (clusterPos[s].z - rhThreadPos.z); float d2 = dist2 / pfClusParams.showerSigma2(); - float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5 * d2); + float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5f * d2); if (rhFracSum[tid] > pfClusParams.minFracTot()) { float fracpct = fraction / rhFracSum[tid]; @@ -531,6 +545,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template >> ALPAKA_FN_ACC static void hcalFastCluster_exotic(const TAcc& acc, reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, + const reco::PFRecHitHCALTopologyDeviceCollection::ConstView topology, int topoId, int nSeeds, int nRHTopo, @@ -572,12 +587,19 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { notDone = true; int i = pfClusteringVars[topoSeedBegin].topoSeedList(); - if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) - rhENormInv = pfClusParams.recHitEnergyNormInvEB_vec()[pfRecHits[i].depth() - 1]; - else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) - rhENormInv = pfClusParams.recHitEnergyNormInvEE_vec()[pfRecHits[i].depth() - 1]; - else - printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + + if (topology.cutsFromDB()) { + rhENormInv = (1.f / topology[pfRecHits[i].denseId()].noiseThreshold()); + } else { + if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) + rhENormInv = pfClusParams.recHitEnergyNormInvHB_vec()[pfRecHits[i].depth() - 1]; + else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) + rhENormInv = pfClusParams.recHitEnergyNormInvHE_vec()[pfRecHits[i].depth() - 1]; + else { + rhENormInv = 0.; + printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + } + } } alpaka::syncBlockThreads(acc); // all threads call sync @@ -651,7 +673,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (clusterPos[s].z - rhThreadPos.z) * (clusterPos[s].z - rhThreadPos.z); float d2 = dist2 / pfClusParams.showerSigma2(); - float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5 * d2); + float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5f * d2); rhFracSum[tid] += fraction; } @@ -669,7 +691,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (clusterPos[s].z - rhThreadPos.z) * (clusterPos[s].z - rhThreadPos.z); float d2 = dist2 / pfClusParams.showerSigma2(); - float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5 * d2); + float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5f * d2); if (rhFracSum[tid] > pfClusParams.minFracTot()) { float fracpct = fraction / rhFracSum[tid]; @@ -798,6 +820,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ALPAKA_FN_ACC static void hcalFastCluster_multiSeedIterative( const TAcc& acc, reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, + const reco::PFRecHitHCALTopologyDeviceCollection::ConstView topology, int topoId, int nSeeds, int nRHTopo, @@ -831,12 +854,19 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { notDone = true; int i = pfClusteringVars[topoSeedBegin].topoSeedList(); - if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) - rhENormInv = pfClusParams.recHitEnergyNormInvEB_vec()[pfRecHits[i].depth() - 1]; - else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) - rhENormInv = pfClusParams.recHitEnergyNormInvEE_vec()[pfRecHits[i].depth() - 1]; - else - printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + + if (topology.cutsFromDB()) { + rhENormInv = (1.f / topology[pfRecHits[i].denseId()].noiseThreshold()); + } else { + if (pfRecHits[i].layer() == PFLayer::HCAL_BARREL1) + rhENormInv = pfClusParams.recHitEnergyNormInvHB_vec()[pfRecHits[i].depth() - 1]; + else if (pfRecHits[i].layer() == PFLayer::HCAL_ENDCAP) + rhENormInv = pfClusParams.recHitEnergyNormInvHE_vec()[pfRecHits[i].depth() - 1]; + else { + rhENormInv = 0.; + printf("Rechit %d has invalid layer %d!\n", i, pfRecHits[i].layer()); + } + } } alpaka::syncBlockThreads(acc); // all threads call sync @@ -910,7 +940,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (clusterPos[s].z - rhThreadPos.z) * (clusterPos[s].z - rhThreadPos.z); float d2 = dist2 / pfClusParams.showerSigma2(); - float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5 * d2); + float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5f * d2); rhFracSum[tid] += fraction; } @@ -928,7 +958,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (clusterPos[s].z - rhThreadPos.z) * (clusterPos[s].z - rhThreadPos.z); float d2 = dist2 / pfClusParams.showerSigma2(); - float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5 * d2); + float fraction = clusterEnergy[s] * rhENormInv * expf(-0.5f * d2); if (rhFracSum[tid] > pfClusParams.minFracTot()) { float fracpct = fraction / rhFracSum[tid]; @@ -1057,6 +1087,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ALPAKA_FN_ACC void operator()(const TAcc& acc, reco::PFClusteringVarsDeviceCollection::View pfClusteringVars, const reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, + const reco::PFRecHitHCALTopologyDeviceCollection::ConstView topology, const reco::PFRecHitHostCollection::ConstView pfRecHits, reco::PFClusterDeviceCollection::View clusterView, reco::PFRecHitFractionDeviceCollection::View fracView, @@ -1082,15 +1113,28 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { int depthOffset = pfRecHits[i].depth() - 1; float energy = pfRecHits[i].energy(); Position3 pos = Position3{pfRecHits[i].x(), pfRecHits[i].y(), pfRecHits[i].z()}; + float seedThreshold = 9999.; + float topoThreshold = 9999.; + + if (topology.cutsFromDB()) { + seedThreshold = topology[pfRecHits[i].denseId()].seedThreshold(); + topoThreshold = topology[pfRecHits[i].denseId()].noiseThreshold(); + } else { + if (layer == PFLayer::HCAL_BARREL1) { + seedThreshold = pfClusParams.seedEThresholdHB_vec()[depthOffset]; + topoThreshold = pfClusParams.topoEThresholdHB_vec()[depthOffset]; + } else if (layer == PFLayer::HCAL_ENDCAP) { + seedThreshold = pfClusParams.seedEThresholdHE_vec()[depthOffset]; + topoThreshold = pfClusParams.topoEThresholdHE_vec()[depthOffset]; + } + } // cmssdt.cern.ch/lxr/source/DataFormats/ParticleFlowReco/interface/PFRecHit.h#0108 float pt2 = energy * energy * (pos.x * pos.x + pos.y * pos.y) / (pos.x * pos.x + pos.y * pos.y + pos.z * pos.z); // Seeding threshold test - if ((layer == PFLayer::HCAL_BARREL1 && energy > pfClusParams.seedEThresholdEB_vec()[depthOffset] && - pt2 > pfClusParams.seedPt2ThresholdEB()) || - (layer == PFLayer::HCAL_ENDCAP && energy > pfClusParams.seedEThresholdEE_vec()[depthOffset] && - pt2 > pfClusParams.seedPt2ThresholdEE())) { + if ((layer == PFLayer::HCAL_BARREL1 && energy > seedThreshold && pt2 > pfClusParams.seedPt2ThresholdHB()) || + (layer == PFLayer::HCAL_ENDCAP && energy > seedThreshold && pt2 > pfClusParams.seedPt2ThresholdHE())) { pfClusteringVars[i].pfrh_isSeed() = 1; for (int k = 0; k < 4; k++) { // Does this seed candidate have a higher energy than four neighbours if (pfRecHits[i].neighbours()(k) < 0) @@ -1104,8 +1148,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { alpaka::atomicAdd(acc, nSeeds, 1u); } // Topo clustering threshold test - if ((layer == PFLayer::HCAL_ENDCAP && energy > pfClusParams.topoEThresholdEE_vec()[depthOffset]) || - (layer == PFLayer::HCAL_BARREL1 && energy > pfClusParams.topoEThresholdEB_vec()[depthOffset])) { + + if ((layer == PFLayer::HCAL_ENDCAP && energy > topoThreshold) || + (layer == PFLayer::HCAL_BARREL1 && energy > topoThreshold)) { pfClusteringVars[i].pfrh_passTopoThresh() = true; pfClusteringVars[i].pfrh_topoId() = i; } else { @@ -1306,6 +1351,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ALPAKA_FN_ACC void operator()(const TAcc& acc, const reco::PFRecHitHostCollection::ConstView pfRecHits, const reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, + const reco::PFRecHitHCALTopologyDeviceCollection::ConstView topology, reco::PFClusteringVarsDeviceCollection::View pfClusteringVars, reco::PFClusterDeviceCollection::View clusterView, reco::PFRecHitFractionDeviceCollection::View fracView) const { @@ -1339,14 +1385,14 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { if (nSeeds == 1) { // Single seed cluster hcalFastCluster_singleSeed( - acc, pfClusParams, topoId, nRHTopo, pfRecHits, pfClusteringVars, clusterView, fracView); + acc, pfClusParams, topology, topoId, nRHTopo, pfRecHits, pfClusteringVars, clusterView, fracView); } else if (nSeeds <= 100 && nRHTopo - nSeeds < threadsPerBlockForClustering) { hcalFastCluster_multiSeedParallel( - acc, pfClusParams, topoId, nSeeds, nRHTopo, pfRecHits, pfClusteringVars, clusterView, fracView); + acc, pfClusParams, topology, topoId, nSeeds, nRHTopo, pfRecHits, pfClusteringVars, clusterView, fracView); } } else if (nSeeds <= 400 && nRHTopo - nSeeds <= 1500) { hcalFastCluster_multiSeedIterative( - acc, pfClusParams, topoId, nSeeds, nRHTopo, pfRecHits, pfClusteringVars, clusterView, fracView); + acc, pfClusParams, topology, topoId, nSeeds, nRHTopo, pfRecHits, pfClusteringVars, clusterView, fracView); } else { if constexpr (debug) { if (once_per_block(acc)) @@ -1367,6 +1413,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ALPAKA_FN_ACC void operator()(const TAcc& acc, const reco::PFRecHitHostCollection::ConstView pfRecHits, const reco::PFClusterParamsDeviceCollection::ConstView pfClusParams, + const reco::PFRecHitHCALTopologyDeviceCollection::ConstView topology, reco::PFClusteringVarsDeviceCollection::View pfClusteringVars, reco::PFClusterDeviceCollection::View clusterView, reco::PFRecHitFractionDeviceCollection::View fracView, @@ -1385,6 +1432,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { if (nRHTopo > 0 && nSeeds > 400 && nRHTopo - nSeeds > 1500) { hcalFastCluster_exotic(acc, pfClusParams, + topology, topoId, nSeeds, nRHTopo, @@ -1420,6 +1468,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { void PFClusterProducerKernel::execute(Queue& queue, const reco::PFClusterParamsDeviceCollection& params, + const reco::PFRecHitHCALTopologyDeviceCollection& topology, reco::PFClusteringVarsDeviceCollection& pfClusteringVars, reco::PFClusteringEdgeVarsDeviceCollection& pfClusteringEdgeVars, const reco::PFRecHitHostCollection& pfRecHits, @@ -1435,6 +1484,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { SeedingTopoThresh{}, pfClusteringVars.view(), params.view(), + topology.view(), pfRecHits.view(), pfClusters.view(), pfrhFractions.view(), @@ -1489,6 +1539,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { FastCluster{}, pfRecHits.view(), params.view(), + topology.view(), pfClusteringVars.view(), pfClusters.view(), pfrhFractions.view()); @@ -1499,6 +1550,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { FastClusterExotic{}, pfRecHits.view(), params.view(), + topology.view(), pfClusteringVars.view(), pfClusters.view(), pfrhFractions.view(), diff --git a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.h b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.h index 8ebafbd366efe..715d484f3120e 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.h +++ b/RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.h @@ -8,6 +8,8 @@ #include "RecoParticleFlow/PFClusterProducer/interface/alpaka/PFClusterParamsDeviceCollection.h" #include "RecoParticleFlow/PFClusterProducer/interface/alpaka/PFClusteringVarsDeviceCollection.h" #include "RecoParticleFlow/PFClusterProducer/interface/alpaka/PFClusteringEdgeVarsDeviceCollection.h" +#include "RecoParticleFlow/PFRecHitProducer/interface/PFRecHitTopologyRecord.h" +#include "RecoParticleFlow/PFRecHitProducer/interface/alpaka/PFRecHitTopologyDeviceCollection.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" namespace ALPAKA_ACCELERATOR_NAMESPACE { @@ -40,6 +42,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { void execute(Queue& queue, const reco::PFClusterParamsDeviceCollection& params, + const reco::PFRecHitHCALTopologyDeviceCollection& topology, reco::PFClusteringVarsDeviceCollection& pfClusteringVars, reco::PFClusteringEdgeVarsDeviceCollection& pfClusteringEdgeVars, const reco::PFRecHitHostCollection& pfRecHits, diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py index 6971ce8a0052a..3436581964004 100644 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py @@ -91,6 +91,6 @@ from Configuration.ProcessModifiers.alpaka_cff import alpaka def _addProcessPFClusterAlpaka(process): - process.load("RecoParticleFlow.PFClusterProducerAlpaka.pfClusterHBHEAlpaka_cff") + process.load("RecoParticleFlow.PFClusterProducer.pfClusterHBHEAlpaka_cff") modifyConfigurationPFClusterAlpaka_ = alpaka.makeProcessModifier(_addProcessPFClusterAlpaka) diff --git a/RecoParticleFlow/PFClusterProducer/python/pfClusterHBHEAlpaka_cff.py b/RecoParticleFlow/PFClusterProducer/python/pfClusterHBHEAlpaka_cff.py index 8053333f0769d..631eee2cec974 100644 --- a/RecoParticleFlow/PFClusterProducer/python/pfClusterHBHEAlpaka_cff.py +++ b/RecoParticleFlow/PFClusterProducer/python/pfClusterHBHEAlpaka_cff.py @@ -61,13 +61,14 @@ pfClusterParamsESProducer = _pfClusterParamsESProducer.clone() pfClusterSoAProducer = _pfClusterSoAProducer.clone( pfRecHits = 'pfRecHitSoAProducerHCAL', + topology = "pfRecHitHCALTopologyESProducer:", pfClusterParams = 'pfClusterParamsESProducer:', synchronise = cms.bool(False) ) legacyPFClusterProducer = _legacyPFClusterProducer.clone( - src = 'pfClusterProducerAlpaka', + src = 'pfClusterSoAProducer', pfClusterParams = 'pfClusterParamsESProducer:', pfClusterBuilder = particleFlowClusterHBHE.pfClusterBuilder, recHitsSource = 'legacyPFRecHitProducer', diff --git a/RecoParticleFlow/PFClusterProducer/test/test_PFRecHitAndClusterSoA.py b/RecoParticleFlow/PFClusterProducer/test/test_PFRecHitAndClusterSoA.py index 4b7e7753cad4e..7cf551884a504 100644 --- a/RecoParticleFlow/PFClusterProducer/test/test_PFRecHitAndClusterSoA.py +++ b/RecoParticleFlow/PFClusterProducer/test/test_PFRecHitAndClusterSoA.py @@ -184,7 +184,7 @@ name = cms.string('PFHBHERecHitCreator'), qualityTests = cms.VPSet( cms.PSet( - usePFThresholdsFromDB = cms.bool(False), + usePFThresholdsFromDB = cms.bool(True), cuts = cms.VPSet( cms.PSet( depth = cms.vint32(1, 2, 3, 4), @@ -224,7 +224,7 @@ name = cms.string('PFHBHERecHitCreator'), qualityTests = cms.VPSet( cms.PSet( - usePFThresholdsFromDB = cms.bool(False), + usePFThresholdsFromDB = cms.bool(True), cuts = cms.VPSet( cms.PSet( depth = cms.vint32(1, 2, 3, 4), @@ -316,7 +316,8 @@ iovIsRunNotTime = cms.bool(True), firstValid = cms.vuint32(1) ) -process.hltParticleFlowRecHitTopologyESProducer = cms.ESProducer(alpaka_backend_str % f"PFRecHit{CAL}TopologyESProducer") +process.hltParticleFlowRecHitTopologyESProducer = cms.ESProducer(alpaka_backend_str % f"PFRecHit{CAL}TopologyESProducer", + usePFThresholdsFromDB = cms.bool(True)) if hcal: process.hltParticleFlowRecHitParamsESProducer = cms.ESProducer(alpaka_backend_str % "PFRecHitHCALParamsESProducer", energyThresholdsHB = cms.vdouble( 0.4, 0.3, 0.3, 0.3 ), @@ -430,6 +431,7 @@ process.hltParticleFlowPFClusterAlpaka = cms.EDProducer(alpaka_backend_str % "PFClusterSoAProducer", pfClusterParams = cms.ESInputTag("hltParticleFlowClusterParamsESProducer:"), + topology = cms.ESInputTag("hltParticleFlowRecHitTopologyESProducer:"), synchronise = cms.bool(args.synchronise)) process.hltParticleFlowPFClusterAlpaka.pfRecHits = cms.InputTag("hltParticleFlowPFRecHitAlpaka") @@ -439,11 +441,12 @@ src = cms.InputTag("hltParticleFlowPFClusterAlpaka"), pfClusterParams = cms.ESInputTag("hltParticleFlowClusterParamsESProducer:"), pfClusterBuilder = process.hltParticleFlowClusterHBHE.pfClusterBuilder, + usePFThresholdsFromDB = cms.bool(True), recHitsSource = cms.InputTag("hltParticleFlowAlpakaToLegacyPFRecHits")) process.hltParticleFlowAlpakaToLegacyPFClusters.PFRecHitsLabelIn = cms.InputTag("hltParticleFlowPFRecHitAlpaka") process.hltParticleFlowClusterHBHE.pfClusterBuilder.maxIterations = 5 -process.hltParticleFlowClusterHBHE.usePFThresholdsFromDB = cms.bool(False) +process.hltParticleFlowClusterHBHE.usePFThresholdsFromDB = cms.bool(True) # Additional customization process.FEVTDEBUGHLToutput.outputCommands = cms.untracked.vstring('drop *_*_*_*') diff --git a/RecoParticleFlow/PFRecHitProducer/interface/PFRecHitTopologyRecord.h b/RecoParticleFlow/PFRecHitProducer/interface/PFRecHitTopologyRecord.h index 78e731ac9eaa6..d539b47ca9ce3 100644 --- a/RecoParticleFlow/PFRecHitProducer/interface/PFRecHitTopologyRecord.h +++ b/RecoParticleFlow/PFRecHitProducer/interface/PFRecHitTopologyRecord.h @@ -5,10 +5,11 @@ #include "FWCore/Framework/interface/EventSetupRecordImplementation.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "Geometry/Records/interface/HcalRecNumberingRecord.h" +#include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h" class PFRecHitHCALTopologyRecord : public edm::eventsetup::DependentRecordImplementation< PFRecHitHCALTopologyRecord, - edm::mpl::Vector> {}; + edm::mpl::Vector> {}; class PFRecHitECALTopologyRecord : public edm::eventsetup::DependentRecordImplementation + + diff --git a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.dev.cc b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.dev.cc index b9287a1f45787..ef18ebc5ecc93 100644 --- a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.dev.cc +++ b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.dev.cc @@ -16,6 +16,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template >> ALPAKA_FN_ACC void operator()(const TAcc& acc, const typename CAL::ParameterType::ConstView params, + const typename CAL::TopologyTypeDevice::ConstView topology, const typename CAL::CaloRecHitSoATypeDevice::ConstView recHits, reco::PFRecHitDeviceCollection::View pfRecHits, uint32_t* __restrict__ denseId2pfRecHit, @@ -23,7 +24,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { // Strided loop over CaloRecHits for (int32_t i : cms::alpakatools::elements_with_stride(acc, recHits.metadata().size())) { // Check energy thresholds/quality cuts (specialised for HCAL/ECAL) - if (!applyCuts(recHits[i], params)) + if (!applyCuts(recHits[i], params, topology)) continue; // Use atomic operation to determine index of the PFRecHit to be constructed @@ -40,7 +41,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } ALPAKA_FN_ACC static bool applyCuts(const typename CAL::CaloRecHitSoATypeDevice::ConstView::const_element rh, - const typename CAL::ParameterType::ConstView params); + const typename CAL::ParameterType::ConstView params, + const typename CAL::TopologyTypeDevice::ConstView topology); ALPAKA_FN_ACC static void constructPFRecHit( reco::PFRecHitDeviceCollection::View::element pfrh, @@ -50,26 +52,33 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template <> ALPAKA_FN_ACC bool PFRecHitProducerKernelConstruct::applyCuts( const typename HCAL::CaloRecHitSoATypeDevice::ConstView::const_element rh, - const HCAL::ParameterType::ConstView params) { + const HCAL::ParameterType::ConstView params, + const HCAL::TopologyTypeDevice::ConstView topology) { // Reject HCAL recHits below enery threshold float threshold = 9999.; const uint32_t detId = rh.detId(); const uint32_t depth = HCAL::getDepth(detId); const uint32_t subdet = getSubdet(detId); - if (subdet == HcalBarrel) { - threshold = params.energyThresholds()[depth - 1]; - } else if (subdet == HcalEndcap) { - threshold = params.energyThresholds()[depth - 1 + HCAL::kMaxDepthHB]; + if (topology.cutsFromDB()) { + threshold = topology.noiseThreshold()[HCAL::detId2denseId(detId)]; } else { - printf("Rechit with detId %u has invalid subdetector %u!\n", detId, subdet); - return false; + if (subdet == HcalBarrel) { + threshold = params.energyThresholds()[depth - 1]; + } else if (subdet == HcalEndcap) { + threshold = params.energyThresholds()[depth - 1 + HCAL::kMaxDepthHB]; + } else { + printf("Rechit with detId %u has invalid subdetector %u!\n", detId, subdet); + return false; + } } return rh.energy() >= threshold; } template <> ALPAKA_FN_ACC bool PFRecHitProducerKernelConstruct::applyCuts( - const ECAL::CaloRecHitSoATypeDevice::ConstView::const_element rh, const ECAL::ParameterType::ConstView params) { + const ECAL::CaloRecHitSoATypeDevice::ConstView::const_element rh, + const ECAL::ParameterType::ConstView params, + const ECAL::TopologyTypeDevice::ConstView topology) { // Reject ECAL recHits below energy threshold if (rh.energy() < params.energyThresholds()[ECAL::detId2denseId(rh.detId())]) return false; @@ -88,6 +97,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { reco::PFRecHitDeviceCollection::View::element pfrh, const HCAL::CaloRecHitSoATypeDevice::ConstView::const_element rh) { pfrh.detId() = rh.detId(); + pfrh.denseId() = HCAL::detId2denseId(rh.detId()); pfrh.energy() = rh.energy(); pfrh.time() = rh.time(); pfrh.depth() = HCAL::getDepth(pfrh.detId()); @@ -105,6 +115,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { reco::PFRecHitDeviceCollection::View::element pfrh, const ECAL::CaloRecHitSoATypeDevice::ConstView::const_element rh) { pfrh.detId() = rh.detId(); + pfrh.denseId() = ECAL::detId2denseId(rh.detId()); pfrh.energy() = rh.energy(); pfrh.time() = rh.time(); pfrh.depth() = 1; @@ -168,11 +179,13 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { void PFRecHitProducerKernel::processRecHits(Queue& queue, const typename CAL::CaloRecHitSoATypeDevice& recHits, const typename CAL::ParameterType& params, + const typename CAL::TopologyTypeDevice& topology, reco::PFRecHitDeviceCollection& pfRecHits) { alpaka::exec(queue, work_div_, PFRecHitProducerKernelConstruct{}, params.view(), + topology.view(), recHits.view(), pfRecHits.view(), denseId2pfRecHit_.data(), diff --git a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.h b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.h index 37638a2370060..ffaef6b0ad748 100644 --- a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.h +++ b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitProducerKernel.h @@ -17,6 +17,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { void processRecHits(Queue& queue, const typename CAL::CaloRecHitSoATypeDevice& recHits, const typename CAL::ParameterType& params, + const typename CAL::TopologyTypeDevice& topology, reco::PFRecHitDeviceCollection& pfRecHits); // Run kernel: Associate topology information (position, neighbours) diff --git a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitSoAProducer.cc b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitSoAProducer.cc index a53ce4f23eed4..8297b6359eaf5 100644 --- a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitSoAProducer.cc +++ b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitSoAProducer.cc @@ -41,7 +41,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { PFRecHitProducerKernel kernel{event.queue(), num_recHits}; for (const auto& token : recHitsToken_) - kernel.processRecHits(event.queue(), event.get(token.first), setup.getData(token.second), pfRecHits); + kernel.processRecHits(event.queue(), event.get(token.first), setup.getData(token.second), topology, pfRecHits); kernel.associateTopologyInfo(event.queue(), topology, pfRecHits); if (synchronise_) diff --git a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitTopologyESProducer.cc b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitTopologyESProducer.cc index c751202f45347..f94db2aecc362 100644 --- a/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitTopologyESProducer.cc +++ b/RecoParticleFlow/PFRecHitProducer/plugins/alpaka/PFRecHitTopologyESProducer.cc @@ -5,6 +5,8 @@ #include #include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h" +#include "CondTools/Hcal/interface/HcalPFCutsHandler.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" @@ -24,15 +26,22 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template class PFRecHitTopologyESProducer : public ESProducer { public: - PFRecHitTopologyESProducer(edm::ParameterSet const& iConfig) : ESProducer(iConfig) { + PFRecHitTopologyESProducer(edm::ParameterSet const& iConfig) + : ESProducer(iConfig), cutsFromDB_(iConfig.getParameter("usePFThresholdsFromDB")) { auto cc = setWhatProduced(this); geomToken_ = cc.consumes(); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { hcalToken_ = cc.consumes(); + hcalCutsToken_ = cc.consumes(); + } } static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; + if constexpr (std::is_same_v) + desc.add("usePFThresholdsFromDB", true); + else // only needs to be true for HBHE + desc.add("usePFThresholdsFromDB", false); descriptions.addWithDefaultLabel(desc); } @@ -40,7 +49,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { const auto& geom = iRecord.get(geomToken_); auto product = std::make_unique(CAL::kSize, cms::alpakatools::host()); auto view = product->view(); - const int calEnums[2] = {CAL::kSubdetectorBarrelId, CAL::kSubdetectorEndcapId}; for (const auto subdet : calEnums) { // Construct topology @@ -61,6 +69,20 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { const uint32_t denseId = CAL::detId2denseId(detId); assert(denseId < CAL::kSize); + // Fill SoA members with HCAL PF Thresholds from GT + if constexpr (std::is_same_v) { + view.cutsFromDB() = false; + if (cutsFromDB_) { + view.cutsFromDB() = true; + const HcalPFCuts& pfCuts = iRecord.get(hcalCutsToken_); + const HcalTopology& htopo = iRecord.get(hcalToken_); + std::unique_ptr prod = std::make_unique(pfCuts); + prod->setTopo(&htopo); + view.noiseThreshold(denseId) = prod->getValues(detId.rawId())->noiseThreshold(); + view.seedThreshold(denseId) = prod->getValues(detId.rawId())->seedThreshold(); + } + } + const GlobalPoint pos = geo->getGeometry(detId)->getPosition(); view.positionX(denseId) = pos.x(); view.positionY(denseId) = pos.y(); @@ -119,6 +141,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { private: edm::ESGetToken geomToken_; edm::ESGetToken hcalToken_; + edm::ESGetToken hcalCutsToken_; + const bool cutsFromDB_; // specialised for HCAL/ECAL, because non-nearest neighbours are defined differently uint32_t getNeighbourDetId(const uint32_t detId, const uint32_t direction, const CaloSubdetectorTopology& topo);