Back to home page

Project CMSSW displayed by LXR

 
 

    


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       // seedFinder
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       // initialClusteringStep
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       // pfClusterBuilder
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 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE
0334 
0335 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h"
0336 DEFINE_FWK_ALPAKA_MODULE(PFClusterSoAProducer);