File indexing completed on 2025-02-05 03:15:15
0001 #include <Eigen/Core>
0002 #include <Eigen/Dense>
0003
0004 #include "DataFormats/ParticleFlowReco/interface/alpaka/PFRecHitDeviceCollection.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/EDGetToken.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010 #include "FWCore/Utilities/interface/StreamID.h"
0011 #include "HeterogeneousCore/AlpakaCore/interface/MoveToDeviceCache.h"
0012 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/SynchronizingEDProducer.h"
0013 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0014 #include "RecoParticleFlow/PFClusterProducer/plugins/alpaka/PFClusterSoAProducerKernel.h"
0015 #include "RecoParticleFlow/PFRecHitProducer/interface/PFRecHitTopologyRecord.h"
0016
0017 #include "PFClusterParamsSoA.h"
0018
0019 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0020 namespace {
0021 using PFClusterParamsCache =
0022 cms::alpakatools::MoveToDeviceCache<Device, PortableHostCollection<::reco::PFClusterParamsSoA>>;
0023 }
0024
0025 class PFClusterSoAProducer : public stream::SynchronizingEDProducer<edm::GlobalCache<PFClusterParamsCache>> {
0026 public:
0027 static std::unique_ptr<PFClusterParamsCache> initializeGlobalCache(edm::ParameterSet const& config) {
0028 constexpr static uint32_t kMaxDepth_barrel = 4;
0029 constexpr static uint32_t kMaxDepth_endcap = 7;
0030 PortableHostCollection<::reco::PFClusterParamsSoA> obj(std::max(kMaxDepth_barrel, kMaxDepth_endcap),
0031 cms::alpakatools::host());
0032 auto view = obj.view();
0033
0034
0035 auto const& sfConf = config.getParameterSet("seedFinder");
0036 view.nNeigh() = sfConf.getParameter<int>("nNeighbours");
0037 auto const& seedFinderConfs = sfConf.getParameterSetVector("thresholdsByDetector");
0038 for (auto const& pset : seedFinderConfs) {
0039 auto const& det = pset.getParameter<std::string>("detector");
0040 auto seedPt2Threshold = std::pow(pset.getParameter<double>("seedingThresholdPt"), 2.);
0041 auto const& thresholds = pset.getParameter<std::vector<double>>("seedingThreshold");
0042 if (det == "HCAL_BARREL1") {
0043 if (thresholds.size() != kMaxDepth_barrel)
0044 throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_barrel
0045 << ") for \"\" vector of det = \"" << det << "\"";
0046 view.seedPt2ThresholdHB() = seedPt2Threshold;
0047 for (size_t idx = 0; idx < thresholds.size(); ++idx) {
0048 view.seedEThresholdHB_vec()[idx] = thresholds[idx];
0049 }
0050 } else if (det == "HCAL_ENDCAP") {
0051 if (thresholds.size() != kMaxDepth_endcap)
0052 throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_endcap
0053 << ") for \"\" vector of det = \"" << det << "\"";
0054 view.seedPt2ThresholdHE() = seedPt2Threshold;
0055 for (size_t idx = 0; idx < thresholds.size(); ++idx) {
0056 view.seedEThresholdHE_vec()[idx] = thresholds[idx];
0057 }
0058 } else {
0059 throw cms::Exception("Configuration") << "Unknown detector when parsing seedFinder: " << det;
0060 }
0061 }
0062
0063
0064 auto const& initConf = config.getParameterSet("initialClusteringStep");
0065 auto const& topoThresholdConf = initConf.getParameterSetVector("thresholdsByDetector");
0066 for (auto const& pset : topoThresholdConf) {
0067 auto const& det = pset.getParameter<std::string>("detector");
0068 auto const& thresholds = pset.getParameter<std::vector<double>>("gatheringThreshold");
0069 if (det == "HCAL_BARREL1") {
0070 if (thresholds.size() != kMaxDepth_barrel)
0071 throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_barrel
0072 << ") for \"\" vector of det = \"" << det << "\"";
0073 for (size_t idx = 0; idx < thresholds.size(); ++idx) {
0074 view.topoEThresholdHB_vec()[idx] = thresholds[idx];
0075 }
0076 } else if (det == "HCAL_ENDCAP") {
0077 if (thresholds.size() != kMaxDepth_endcap)
0078 throw cms::Exception("Configuration") << "Invalid size (" << thresholds.size() << " != " << kMaxDepth_endcap
0079 << ") for \"\" vector of det = \"" << det << "\"";
0080 for (size_t idx = 0; idx < thresholds.size(); ++idx) {
0081 view.topoEThresholdHE_vec()[idx] = thresholds[idx];
0082 }
0083 } else {
0084 throw cms::Exception("Configuration") << "Unknown detector when parsing initClusteringStep: " << det;
0085 }
0086 }
0087
0088
0089 auto const& pfClusterPSet = config.getParameterSet("pfClusterBuilder");
0090 view.showerSigma2() = std::pow(pfClusterPSet.getParameter<double>("showerSigma"), 2.);
0091 view.minFracToKeep() = pfClusterPSet.getParameter<double>("minFractionToKeep");
0092 view.minFracTot() = pfClusterPSet.getParameter<double>("minFracTot");
0093 view.maxIterations() = pfClusterPSet.getParameter<unsigned int>("maxIterations");
0094 view.excludeOtherSeeds() = pfClusterPSet.getParameter<bool>("excludeOtherSeeds");
0095 view.stoppingTolerance() = pfClusterPSet.getParameter<double>("stoppingTolerance");
0096 auto const& pcPSet = pfClusterPSet.getParameterSet("positionCalc");
0097 view.minFracInCalc() = pcPSet.getParameter<double>("minFractionInCalc");
0098 view.minAllowedNormalization() = pcPSet.getParameter<double>("minAllowedNormalization");
0099
0100 auto const& recHitEnergyNormConf = pfClusterPSet.getParameterSetVector("recHitEnergyNorms");
0101 for (auto const& pset : recHitEnergyNormConf) {
0102 auto const& recHitNorms = pset.getParameter<std::vector<double>>("recHitEnergyNorm");
0103 auto const& det = pset.getParameter<std::string>("detector");
0104 if (det == "HCAL_BARREL1") {
0105 if (recHitNorms.size() != kMaxDepth_barrel)
0106 throw cms::Exception("Configuration")
0107 << "Invalid size (" << recHitNorms.size() << " != " << kMaxDepth_barrel
0108 << ") for \"\" vector of det = \"" << det << "\"";
0109 for (size_t idx = 0; idx < recHitNorms.size(); ++idx) {
0110 view.recHitEnergyNormInvHB_vec()[idx] = 1. / recHitNorms[idx];
0111 }
0112 } else if (det == "HCAL_ENDCAP") {
0113 if (recHitNorms.size() != kMaxDepth_endcap)
0114 throw cms::Exception("Configuration")
0115 << "Invalid size (" << recHitNorms.size() << " != " << kMaxDepth_endcap
0116 << ") for \"\" vector of det = \"" << det << "\"";
0117 for (size_t idx = 0; idx < recHitNorms.size(); ++idx) {
0118 view.recHitEnergyNormInvHE_vec()[idx] = 1. / recHitNorms[idx];
0119 }
0120 } else {
0121 throw cms::Exception("Configuration") << "Unknown detector when parsing recHitEnergyNorms: " << det;
0122 }
0123 }
0124
0125 auto const& barrelTimeResConf = pfClusterPSet.getParameterSet("timeResolutionCalcBarrel");
0126 view.barrelTimeResConsts_corrTermLowE() = barrelTimeResConf.getParameter<double>("corrTermLowE");
0127 view.barrelTimeResConsts_threshLowE() = barrelTimeResConf.getParameter<double>("threshLowE");
0128 view.barrelTimeResConsts_noiseTerm() = barrelTimeResConf.getParameter<double>("noiseTerm");
0129 view.barrelTimeResConsts_constantTermLowE2() =
0130 std::pow(barrelTimeResConf.getParameter<double>("constantTermLowE"), 2.);
0131 view.barrelTimeResConsts_noiseTermLowE() = barrelTimeResConf.getParameter<double>("noiseTermLowE");
0132 view.barrelTimeResConsts_threshHighE() = barrelTimeResConf.getParameter<double>("threshHighE");
0133 view.barrelTimeResConsts_constantTerm2() = std::pow(barrelTimeResConf.getParameter<double>("constantTerm"), 2.);
0134 view.barrelTimeResConsts_resHighE2() =
0135 std::pow(view.barrelTimeResConsts_noiseTerm() / view.barrelTimeResConsts_threshHighE(), 2.) +
0136 view.barrelTimeResConsts_constantTerm2();
0137
0138 auto const& endcapTimeResConf = pfClusterPSet.getParameterSet("timeResolutionCalcEndcap");
0139 view.endcapTimeResConsts_corrTermLowE() = endcapTimeResConf.getParameter<double>("corrTermLowE");
0140 view.endcapTimeResConsts_threshLowE() = endcapTimeResConf.getParameter<double>("threshLowE");
0141 view.endcapTimeResConsts_noiseTerm() = endcapTimeResConf.getParameter<double>("noiseTerm");
0142 view.endcapTimeResConsts_constantTermLowE2() =
0143 std::pow(endcapTimeResConf.getParameter<double>("constantTermLowE"), 2.);
0144 view.endcapTimeResConsts_noiseTermLowE() = endcapTimeResConf.getParameter<double>("noiseTermLowE");
0145 view.endcapTimeResConsts_threshHighE() = endcapTimeResConf.getParameter<double>("threshHighE");
0146 view.endcapTimeResConsts_constantTerm2() = std::pow(endcapTimeResConf.getParameter<double>("constantTerm"), 2.);
0147 view.endcapTimeResConsts_resHighE2() =
0148 std::pow(view.endcapTimeResConsts_noiseTerm() / view.endcapTimeResConsts_threshHighE(), 2.) +
0149 view.endcapTimeResConsts_constantTerm2();
0150
0151 return std::make_unique<PFClusterParamsCache>(std::move(obj));
0152 }
0153
0154 PFClusterSoAProducer(edm::ParameterSet const& config, PFClusterParamsCache const*)
0155 : SynchronizingEDProducer(config),
0156 topologyToken_(esConsumes(config.getParameter<edm::ESInputTag>("topology"))),
0157 inputPFRecHitSoA_Token_{consumes(config.getParameter<edm::InputTag>("pfRecHits"))},
0158 inputPFRecHitNum_Token_{consumes(config.getParameter<edm::InputTag>("pfRecHits"))},
0159 outputPFClusterSoA_Token_{produces()},
0160 outputPFRHFractionSoA_Token_{produces()},
0161 numRHF_{cms::alpakatools::make_host_buffer<uint32_t, Platform>()},
0162 synchronise_(config.getParameter<bool>("synchronise")) {}
0163
0164 void acquire(device::Event const& event, device::EventSetup const& setup) override {
0165 const reco::PFRecHitHCALTopologyDeviceCollection& topology = setup.getData(topologyToken_);
0166 const reco::PFRecHitDeviceCollection& pfRecHits = event.get(inputPFRecHitSoA_Token_);
0167 int nRH = event.get(inputPFRecHitNum_Token_);
0168 const auto& params = globalCache()->get(event.queue());
0169
0170 pfClusteringVars_.emplace(nRH, event.queue());
0171 pfClusteringEdgeVars_.emplace(nRH * 8, event.queue());
0172 pfClusters_.emplace(nRH, event.queue());
0173
0174 *numRHF_ = 0;
0175
0176 if (nRH != 0) {
0177 PFClusterProducerKernel kernel(event.queue());
0178 kernel.seedTopoAndContract(event.queue(),
0179 params.const_view(),
0180 topology,
0181 *pfClusteringVars_,
0182 *pfClusteringEdgeVars_,
0183 pfRecHits,
0184 nRH,
0185 *pfClusters_,
0186 numRHF_.data());
0187 }
0188 }
0189
0190 void produce(device::Event& event, device::EventSetup const& setup) override {
0191 const auto& params = globalCache()->get(event.queue());
0192 const reco::PFRecHitHCALTopologyDeviceCollection& topology = setup.getData(topologyToken_);
0193 const reco::PFRecHitDeviceCollection& pfRecHits = event.get(inputPFRecHitSoA_Token_);
0194
0195 std::optional<reco::PFRecHitFractionDeviceCollection> pfrhFractions;
0196
0197 int nRH = event.get(inputPFRecHitNum_Token_);
0198 if (nRH != 0) {
0199 pfrhFractions.emplace(*numRHF_, event.queue());
0200 PFClusterProducerKernel kernel(event.queue());
0201 kernel.cluster(event.queue(),
0202 params.const_view(),
0203 topology,
0204 *pfClusteringVars_,
0205 *pfClusteringEdgeVars_,
0206 pfRecHits,
0207 nRH,
0208 *pfClusters_,
0209 *pfrhFractions);
0210 } else {
0211 pfrhFractions.emplace(0, event.queue());
0212 }
0213
0214 if (synchronise_)
0215 alpaka::wait(event.queue());
0216
0217 event.emplace(outputPFClusterSoA_Token_, std::move(*pfClusters_));
0218 event.emplace(outputPFRHFractionSoA_Token_, std::move(*pfrhFractions));
0219 }
0220
0221 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0222 edm::ParameterSetDescription desc;
0223 desc.add<edm::InputTag>("pfRecHits", edm::InputTag(""));
0224 desc.add<edm::ESInputTag>("topology", edm::ESInputTag(""));
0225
0226 {
0227 auto const psetName = "seedFinder";
0228 edm::ParameterSetDescription foo;
0229 foo.add<int>("nNeighbours", 4);
0230 {
0231 edm::ParameterSetDescription validator;
0232 validator.add<std::string>("detector", "");
0233 validator.add<std::vector<double>>("seedingThreshold", {});
0234 validator.add<double>("seedingThresholdPt", 0.);
0235 std::vector<edm::ParameterSet> vDefaults(2);
0236 vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0237 vDefaults[0].addParameter<std::vector<double>>("seedingThreshold", {0.125, 0.25, 0.35, 0.35});
0238 vDefaults[0].addParameter<double>("seedingThresholdPt", 0.);
0239 vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0240 vDefaults[1].addParameter<std::vector<double>>("seedingThreshold",
0241 {0.1375, 0.275, 0.275, 0.275, 0.275, 0.275, 0.275});
0242 vDefaults[1].addParameter<double>("seedingThresholdPt", 0.);
0243 foo.addVPSet("thresholdsByDetector", validator, vDefaults);
0244 }
0245 desc.add(psetName, foo);
0246 }
0247 {
0248 auto const psetName = "initialClusteringStep";
0249 edm::ParameterSetDescription foo;
0250 {
0251 edm::ParameterSetDescription validator;
0252 validator.add<std::string>("detector", "");
0253 validator.add<std::vector<double>>("gatheringThreshold", {});
0254 std::vector<edm::ParameterSet> vDefaults(2);
0255 vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0256 vDefaults[0].addParameter<std::vector<double>>("gatheringThreshold", {0.1, 0.2, 0.3, 0.3});
0257 vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0258 vDefaults[1].addParameter<std::vector<double>>("gatheringThreshold", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0259 foo.addVPSet("thresholdsByDetector", validator, vDefaults);
0260 }
0261 desc.add(psetName, foo);
0262 }
0263 {
0264 auto const psetName = "pfClusterBuilder";
0265 edm::ParameterSetDescription foo;
0266 foo.add<unsigned int>("maxIterations", 50);
0267 foo.add<double>("minFracTot", 1e-20);
0268 foo.add<double>("minFractionToKeep", 1e-7);
0269 foo.add<bool>("excludeOtherSeeds", true);
0270 foo.add<double>("showerSigma", 10.);
0271 foo.add<double>("stoppingTolerance", 1e-8);
0272 {
0273 edm::ParameterSetDescription validator;
0274 validator.add<std::string>("detector", "");
0275 validator.add<std::vector<double>>("recHitEnergyNorm", {});
0276 std::vector<edm::ParameterSet> vDefaults(2);
0277 vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0278 vDefaults[0].addParameter<std::vector<double>>("recHitEnergyNorm", {0.1, 0.2, 0.3, 0.3});
0279 vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0280 vDefaults[1].addParameter<std::vector<double>>("recHitEnergyNorm", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0281 foo.addVPSet("recHitEnergyNorms", validator, vDefaults);
0282 }
0283 {
0284 edm::ParameterSetDescription bar;
0285 bar.add<double>("minFractionInCalc", 1e-9);
0286 bar.add<double>("minAllowedNormalization", 1e-9);
0287 foo.add("positionCalc", bar);
0288 }
0289 {
0290 edm::ParameterSetDescription bar;
0291 bar.add<double>("corrTermLowE", 0.);
0292 bar.add<double>("threshLowE", 6.);
0293 bar.add<double>("noiseTerm", 21.86);
0294 bar.add<double>("constantTermLowE", 4.24);
0295 bar.add<double>("noiseTermLowE", 8.);
0296 bar.add<double>("threshHighE", 15.);
0297 bar.add<double>("constantTerm", 2.82);
0298 foo.add("timeResolutionCalcBarrel", bar);
0299 }
0300 {
0301 edm::ParameterSetDescription bar;
0302 bar.add<double>("corrTermLowE", 0.);
0303 bar.add<double>("threshLowE", 6.);
0304 bar.add<double>("noiseTerm", 21.86);
0305 bar.add<double>("constantTermLowE", 4.24);
0306 bar.add<double>("noiseTermLowE", 8.);
0307 bar.add<double>("threshHighE", 15.);
0308 bar.add<double>("constantTerm", 2.82);
0309 foo.add("timeResolutionCalcEndcap", bar);
0310 }
0311 desc.add(psetName, foo);
0312 }
0313
0314 desc.add<bool>("synchronise", false);
0315 descriptions.addWithDefaultLabel(desc);
0316 }
0317
0318 static void globalEndJob(PFClusterParamsCache*) {}
0319
0320 private:
0321 const device::ESGetToken<reco::PFRecHitHCALTopologyDeviceCollection, PFRecHitHCALTopologyRecord> topologyToken_;
0322 const device::EDGetToken<reco::PFRecHitDeviceCollection> inputPFRecHitSoA_Token_;
0323 const edm::EDGetTokenT<cms_uint32_t> inputPFRecHitNum_Token_;
0324 const device::EDPutToken<reco::PFClusterDeviceCollection> outputPFClusterSoA_Token_;
0325 const device::EDPutToken<reco::PFRecHitFractionDeviceCollection> outputPFRHFractionSoA_Token_;
0326 cms::alpakatools::host_buffer<uint32_t> numRHF_;
0327 std::optional<reco::PFClusteringVarsDeviceCollection> pfClusteringVars_;
0328 std::optional<reco::PFClusteringEdgeVarsDeviceCollection> pfClusteringEdgeVars_;
0329 std::optional<reco::PFClusterDeviceCollection> pfClusters_;
0330 const bool synchronise_;
0331 };
0332
0333 }
0334
0335 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h"
0336 DEFINE_FWK_ALPAKA_MODULE(PFClusterSoAProducer);