File indexing completed on 2024-06-22 02:24:04
0001 #include <Eigen/Core>
0002 #include <Eigen/Dense>
0003 #include <algorithm>
0004 #include <array>
0005 #include <cmath>
0006 #include <iostream>
0007 #include <memory>
0008 #include <string>
0009 #include <type_traits>
0010 #include <unordered_map>
0011 #include <utility>
0012 #include <vector>
0013
0014 #include "FWCore/Framework/interface/stream/EDProducer.h"
0015 #include "FWCore/Utilities/interface/StreamID.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/Utilities/interface/EDGetToken.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Utilities/interface/EDPutToken.h"
0024
0025 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0026 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0027
0028 #include "DataFormats/Common/interface/Handle.h"
0029 #include "DataFormats/DetId/interface/DetId.h"
0030 #include "DataFormats/ParticleFlowReco/interface/PFRecHitHostCollection.h"
0031 #include "DataFormats/ParticleFlowReco/interface/PFClusterHostCollection.h"
0032 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFractionHostCollection.h"
0033 #include "HeterogeneousCore/CUDACore/interface/JobConfigurationGPURecord.h"
0034 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0035
0036 class LegacyPFClusterProducer : public edm::stream::EDProducer<> {
0037 public:
0038 LegacyPFClusterProducer(edm::ParameterSet const& config)
0039 : pfClusterSoAToken_(consumes(config.getParameter<edm::InputTag>("src"))),
0040 pfRecHitFractionSoAToken_(consumes(config.getParameter<edm::InputTag>("src"))),
0041 InputPFRecHitSoA_Token_{consumes(config.getParameter<edm::InputTag>("PFRecHitsLabelIn"))},
0042 legacyPfClustersToken_(produces()),
0043 recHitsLabel_(consumes(config.getParameter<edm::InputTag>("recHitsSource"))),
0044 hcalCutsToken_(esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"))),
0045 cutsFromDB_(config.getParameter<bool>("usePFThresholdsFromDB")) {
0046 edm::ConsumesCollector cc = consumesCollector();
0047
0048
0049 const edm::ParameterSet& pfcConf = config.getParameterSet("pfClusterBuilder");
0050 if (!pfcConf.empty()) {
0051 const auto& acConf = pfcConf.getParameterSet("positionCalc");
0052 if (!acConf.empty()) {
0053 const std::string& algoac = acConf.getParameter<std::string>("algoName");
0054 if (!algoac.empty())
0055 positionCalc_ = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0056 }
0057
0058 const auto& acConf2 = pfcConf.getParameterSet("allCellsPositionCalc");
0059 if (!acConf2.empty()) {
0060 const std::string& algoac = acConf2.getParameter<std::string>("algoName");
0061 if (!algoac.empty())
0062 allCellsPositionCalc_ = PFCPositionCalculatorFactory::get()->create(algoac, acConf2, cc);
0063 }
0064 }
0065 }
0066
0067 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0068 edm::ParameterSetDescription desc;
0069 desc.add<edm::InputTag>("src", edm::InputTag("pfClusterSoAProducer"));
0070 desc.add<edm::InputTag>("PFRecHitsLabelIn", edm::InputTag("pfRecHitSoAProducerHCAL"));
0071 desc.add<edm::InputTag>("recHitsSource", edm::InputTag("legacyPFRecHitProducer"));
0072 desc.add<bool>("usePFThresholdsFromDB", true);
0073 {
0074 edm::ParameterSetDescription pfClusterBuilder;
0075 pfClusterBuilder.add<unsigned int>("maxIterations", 5);
0076 pfClusterBuilder.add<double>("minFracTot", 1e-20);
0077 pfClusterBuilder.add<double>("minFractionToKeep", 1e-7);
0078 pfClusterBuilder.add<bool>("excludeOtherSeeds", true);
0079 pfClusterBuilder.add<double>("showerSigma", 10.);
0080 pfClusterBuilder.add<double>("stoppingTolerance", 1e-8);
0081 pfClusterBuilder.add<double>("timeSigmaEB", 10.);
0082 pfClusterBuilder.add<double>("timeSigmaEE", 10.);
0083 pfClusterBuilder.add<double>("maxNSigmaTime", 10.);
0084 pfClusterBuilder.add<double>("minChi2Prob", 0.);
0085 pfClusterBuilder.add<bool>("clusterTimeResFromSeed", false);
0086 pfClusterBuilder.add<std::string>("algoName", "");
0087 pfClusterBuilder.add<edm::ParameterSetDescription>("positionCalcForConvergence", {});
0088 {
0089 edm::ParameterSetDescription validator;
0090 validator.add<std::string>("detector", "");
0091 validator.add<std::vector<int>>("depths", {});
0092 validator.add<std::vector<double>>("recHitEnergyNorm", {});
0093 std::vector<edm::ParameterSet> vDefaults(2);
0094 vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0095 vDefaults[0].addParameter<std::vector<int>>("depths", {1, 2, 3, 4});
0096 vDefaults[0].addParameter<std::vector<double>>("recHitEnergyNorm", {0.1, 0.2, 0.3, 0.3});
0097 vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0098 vDefaults[1].addParameter<std::vector<int>>("depths", {1, 2, 3, 4, 5, 6, 7});
0099 vDefaults[1].addParameter<std::vector<double>>("recHitEnergyNorm", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0100 pfClusterBuilder.addVPSet("recHitEnergyNorms", validator, vDefaults);
0101 }
0102 {
0103 edm::ParameterSetDescription bar;
0104 bar.add<std::string>("algoName", "Basic2DGenericPFlowPositionCalc");
0105 bar.add<double>("minFractionInCalc", 1e-9);
0106 bar.add<int>("posCalcNCrystals", 5);
0107 {
0108 edm::ParameterSetDescription validator;
0109 validator.add<std::string>("detector", "");
0110 validator.add<std::vector<int>>("depths", {});
0111 validator.add<std::vector<double>>("logWeightDenominator", {});
0112 std::vector<edm::ParameterSet> vDefaults(2);
0113 vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0114 vDefaults[0].addParameter<std::vector<int>>("depths", {1, 2, 3, 4});
0115 vDefaults[0].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.3, 0.3});
0116 vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0117 vDefaults[1].addParameter<std::vector<int>>("depths", {1, 2, 3, 4, 5, 6, 7});
0118 vDefaults[1].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0119 bar.addVPSet("logWeightDenominatorByDetector", validator, vDefaults);
0120 }
0121 bar.add<double>("minAllowedNormalization", 1e-9);
0122 bar.add<edm::ParameterSetDescription>("timeResolutionCalcBarrel", {});
0123 bar.add<edm::ParameterSetDescription>("timeResolutionCalcEndcap", {});
0124 pfClusterBuilder.add("positionCalc", bar);
0125 }
0126 {
0127 edm::ParameterSetDescription bar;
0128 bar.add<std::string>("algoName", "Basic2DGenericPFlowPositionCalc");
0129 bar.add<double>("minFractionInCalc", 1e-9);
0130 bar.add<int>("posCalcNCrystals", -1);
0131 {
0132 edm::ParameterSetDescription validator;
0133 validator.add<std::string>("detector", "");
0134 validator.add<std::vector<int>>("depths", {});
0135 validator.add<std::vector<double>>("logWeightDenominator", {});
0136 std::vector<edm::ParameterSet> vDefaults(2);
0137 vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0138 vDefaults[0].addParameter<std::vector<int>>("depths", {1, 2, 3, 4});
0139 vDefaults[0].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.3, 0.3});
0140 vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0141 vDefaults[1].addParameter<std::vector<int>>("depths", {1, 2, 3, 4, 5, 6, 7});
0142 vDefaults[1].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0143 bar.addVPSet("logWeightDenominatorByDetector", validator, vDefaults);
0144 }
0145 bar.add<double>("minAllowedNormalization", 1e-9);
0146 bar.add<edm::ParameterSetDescription>("timeResolutionCalcBarrel", {});
0147 bar.add<edm::ParameterSetDescription>("timeResolutionCalcEndcap", {});
0148 pfClusterBuilder.add("allCellsPositionCalc", bar);
0149 }
0150 {
0151 edm::ParameterSetDescription bar;
0152 bar.add<double>("corrTermLowE", 0.);
0153 bar.add<double>("threshLowE", 6.);
0154 bar.add<double>("noiseTerm", 21.86);
0155 bar.add<double>("constantTermLowE", 4.24);
0156 bar.add<double>("noiseTermLowE", 8.);
0157 bar.add<double>("threshHighE", 15.);
0158 bar.add<double>("constantTerm", 2.82);
0159 pfClusterBuilder.add("timeResolutionCalcBarrel", bar);
0160 }
0161 {
0162 edm::ParameterSetDescription bar;
0163 bar.add<double>("corrTermLowE", 0.);
0164 bar.add<double>("threshLowE", 6.);
0165 bar.add<double>("noiseTerm", 21.86);
0166 bar.add<double>("constantTermLowE", 4.24);
0167 bar.add<double>("noiseTermLowE", 8.);
0168 bar.add<double>("threshHighE", 15.);
0169 bar.add<double>("constantTerm", 2.82);
0170 pfClusterBuilder.add("timeResolutionCalcEndcap", bar);
0171 }
0172 {
0173 edm::ParameterSetDescription bar;
0174 pfClusterBuilder.add("positionReCalc", bar);
0175 }
0176 {
0177 edm::ParameterSetDescription bar;
0178 pfClusterBuilder.add("energyCorrector", bar);
0179 }
0180 desc.add("pfClusterBuilder", pfClusterBuilder);
0181 }
0182 descriptions.addWithDefaultLabel(desc);
0183 }
0184
0185 private:
0186 void produce(edm::Event&, const edm::EventSetup&) override;
0187 const edm::EDGetTokenT<reco::PFClusterHostCollection> pfClusterSoAToken_;
0188 const edm::EDGetTokenT<reco::PFRecHitFractionHostCollection> pfRecHitFractionSoAToken_;
0189 const edm::EDGetTokenT<reco::PFRecHitHostCollection> InputPFRecHitSoA_Token_;
0190 const edm::EDPutTokenT<reco::PFClusterCollection> legacyPfClustersToken_;
0191 const edm::EDGetTokenT<reco::PFRecHitCollection> recHitsLabel_;
0192 const edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0193 const bool cutsFromDB_;
0194
0195 std::unique_ptr<PFCPositionCalculatorBase> positionCalc_;
0196 std::unique_ptr<PFCPositionCalculatorBase> allCellsPositionCalc_;
0197 };
0198
0199 void LegacyPFClusterProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0200 const reco::PFRecHitHostCollection& pfRecHits = event.get(InputPFRecHitSoA_Token_);
0201
0202 HcalPFCuts const* paramPF = cutsFromDB_ ? &setup.getData(hcalCutsToken_) : nullptr;
0203
0204 auto const& pfClusterSoA = event.get(pfClusterSoAToken_).const_view();
0205 auto const& pfRecHitFractionSoA = event.get(pfRecHitFractionSoAToken_).const_view();
0206
0207 int nRH = 0;
0208 if (pfRecHits->metadata().size() != 0)
0209 nRH = pfRecHits.view().size();
0210 reco::PFClusterCollection out;
0211 out.reserve(nRH);
0212
0213 auto const rechitsHandle = event.getHandle(recHitsLabel_);
0214
0215 if (nRH != 0) {
0216
0217 std::vector<int> nTopoSeeds(nRH, 0);
0218
0219 for (int i = 0; i < pfClusterSoA.nSeeds(); i++) {
0220 nTopoSeeds[pfClusterSoA[i].topoId()]++;
0221 }
0222
0223
0224 for (int i = 0; i < pfClusterSoA.nSeeds(); i++) {
0225 unsigned int n = pfClusterSoA[i].seedRHIdx();
0226 reco::PFCluster temp;
0227 temp.setSeed((*rechitsHandle)[n].detId());
0228 int offset = pfClusterSoA[i].rhfracOffset();
0229 for (int k = offset; k < (offset + pfClusterSoA[i].rhfracSize()) && k >= 0;
0230 k++) {
0231 if (pfRecHitFractionSoA[k].pfrhIdx() < nRH && pfRecHitFractionSoA[k].pfrhIdx() > -1 &&
0232 pfRecHitFractionSoA[k].frac() > 0.0) {
0233 const reco::PFRecHitRef& refhit = reco::PFRecHitRef(rechitsHandle, pfRecHitFractionSoA[k].pfrhIdx());
0234 temp.addRecHitFraction(reco::PFRecHitFraction(refhit, pfRecHitFractionSoA[k].frac()));
0235 }
0236 }
0237
0238
0239 if (nTopoSeeds[pfClusterSoA[i].topoId()] == 1 && allCellsPositionCalc_) {
0240 allCellsPositionCalc_->calculateAndSetPosition(temp, paramPF);
0241 } else {
0242 positionCalc_->calculateAndSetPosition(temp, paramPF);
0243 }
0244 out.emplace_back(std::move(temp));
0245 }
0246 }
0247
0248 event.emplace(legacyPfClustersToken_, std::move(out));
0249 }
0250
0251 #include "FWCore/Framework/interface/MakerMacros.h"
0252 DEFINE_FWK_MODULE(LegacyPFClusterProducer);
