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