File indexing completed on 2024-06-22 02:24:04
0001 #include "DataFormats/Math/interface/deltaPhi.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0004 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "Math/GenVector/VectorUtil.h"
0008 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterBuilderBase.h"
0009 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0010 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0011 #include "vdt/vdtMath.h"
0012
0013 #include <iterator>
0014 #include <unordered_map>
0015
0016 class PFMultiDepthClusterizer final : public PFClusterBuilderBase {
0017 typedef PFMultiDepthClusterizer B2DGPF;
0018
0019 public:
0020 PFMultiDepthClusterizer(const edm::ParameterSet& conf, edm::ConsumesCollector& cc);
0021
0022 ~PFMultiDepthClusterizer() override = default;
0023 PFMultiDepthClusterizer(const B2DGPF&) = delete;
0024 B2DGPF& operator=(const B2DGPF&) = delete;
0025
0026 void update(const edm::EventSetup& es) override { _allCellsPosCalc->update(es); }
0027
0028 void buildClusters(const reco::PFClusterCollection&,
0029 const std::vector<bool>&,
0030 reco::PFClusterCollection& outclus,
0031 const HcalPFCuts*) override;
0032
0033 private:
0034 std::unique_ptr<PFCPositionCalculatorBase> _allCellsPosCalc;
0035 double nSigmaEta_;
0036 double nSigmaPhi_;
0037
0038 class ClusterLink {
0039 public:
0040 ClusterLink(unsigned int i, unsigned int j, double DR, int DZ, double energy) {
0041 from_ = i;
0042 to_ = j;
0043 linkDR_ = DR;
0044 linkDZ_ = DZ;
0045 linkE_ = energy;
0046 }
0047
0048 ~ClusterLink() = default;
0049
0050 unsigned int from() const { return from_; }
0051 unsigned int to() const { return to_; }
0052 double dR() const { return linkDR_; }
0053 int dZ() const { return linkDZ_; }
0054 double energy() const { return linkE_; }
0055
0056 private:
0057 unsigned int from_;
0058 unsigned int to_;
0059 double linkDR_;
0060 int linkDZ_;
0061 double linkE_;
0062 };
0063
0064 void calculateShowerShapes(const reco::PFClusterCollection&, std::vector<double>&, std::vector<double>&);
0065 std::vector<ClusterLink> link(const reco::PFClusterCollection&,
0066 const std::vector<double>&,
0067 const std::vector<double>&);
0068 std::vector<ClusterLink> prune(std::vector<ClusterLink>&, std::vector<bool>& linkedClusters);
0069
0070 void expandCluster(reco::PFCluster&,
0071 unsigned int point,
0072 std::vector<bool>& mask,
0073 const reco::PFClusterCollection&,
0074 const std::vector<ClusterLink>& links);
0075
0076 void absorbCluster(reco::PFCluster&, const reco::PFCluster&);
0077 };
0078
0079 DEFINE_EDM_PLUGIN(PFClusterBuilderFactory, PFMultiDepthClusterizer, "PFMultiDepthClusterizer");
0080
0081 PFMultiDepthClusterizer::PFMultiDepthClusterizer(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0082 : PFClusterBuilderBase(conf, cc) {
0083 const auto& acConf = conf.getParameterSet("allCellsPositionCalc");
0084 if (!acConf.empty()) {
0085 const std::string& algoac = acConf.getParameter<std::string>("algoName");
0086 if (!algoac.empty())
0087 _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0088 }
0089
0090 nSigmaEta_ = pow(conf.getParameter<double>("nSigmaEta"), 2);
0091 nSigmaPhi_ = pow(conf.getParameter<double>("nSigmaPhi"), 2);
0092 }
0093
0094 void PFMultiDepthClusterizer::buildClusters(const reco::PFClusterCollection& input,
0095 const std::vector<bool>& seedable,
0096 reco::PFClusterCollection& output,
0097 const HcalPFCuts* hcalCuts) {
0098 std::vector<double> etaRMS2(input.size(), 0.0);
0099 std::vector<double> phiRMS2(input.size(), 0.0);
0100
0101
0102
0103
0104
0105
0106 calculateShowerShapes(input, etaRMS2, phiRMS2);
0107
0108
0109 auto&& links = link(input, etaRMS2, phiRMS2);
0110
0111
0112
0113 std::vector<bool> mask(input.size(), false);
0114 std::vector<bool> linked(input.size(), false);
0115
0116
0117 auto&& prunedLinks = prune(links, linked);
0118
0119
0120
0121
0122
0123
0124 for (unsigned int i = 0; i < input.size(); ++i) {
0125
0126 if (mask[i])
0127 continue;
0128
0129 if (!linked[i]) {
0130 output.push_back(input[i]);
0131
0132 mask[i] = true;
0133 continue;
0134 }
0135
0136
0137 reco::PFCluster cluster = input[i];
0138 mask[i] = true;
0139 expandCluster(cluster, i, mask, input, prunedLinks);
0140 _allCellsPosCalc->calculateAndSetPosition(cluster, hcalCuts);
0141 output.push_back(cluster);
0142
0143 }
0144 }
0145
0146 void PFMultiDepthClusterizer::calculateShowerShapes(const reco::PFClusterCollection& clusters,
0147 std::vector<double>& etaRMS2,
0148 std::vector<double>& phiRMS2) {
0149
0150
0151 for (unsigned int i = 0; i < clusters.size(); ++i) {
0152 const reco::PFCluster& cluster = clusters[i];
0153 double etaSum = 0.0;
0154 double phiSum = 0.0;
0155 auto const& crep = cluster.positionREP();
0156 for (const auto& frac : cluster.recHitFractions()) {
0157 auto const& h = *frac.recHitRef();
0158 auto const& rep = h.positionREP();
0159 etaSum += (frac.fraction() * h.energy()) * std::abs(rep.eta() - crep.eta());
0160 phiSum += (frac.fraction() * h.energy()) * std::abs(deltaPhi(rep.phi(), crep.phi()));
0161 }
0162
0163 etaRMS2[i] = std::max(etaSum / cluster.energy(), 0.1);
0164 etaRMS2[i] *= etaRMS2[i];
0165 phiRMS2[i] = std::max(phiSum / cluster.energy(), 0.1);
0166 phiRMS2[i] *= phiRMS2[i];
0167 }
0168 }
0169
0170 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::link(
0171 const reco::PFClusterCollection& clusters, const std::vector<double>& etaRMS2, const std::vector<double>& phiRMS2) {
0172 std::vector<ClusterLink> links;
0173
0174 for (unsigned int i = 0; i < clusters.size(); ++i)
0175 for (unsigned int j = 0; j < clusters.size(); ++j) {
0176 if (i == j)
0177 continue;
0178
0179 const reco::PFCluster& cluster1 = clusters[i];
0180 const reco::PFCluster& cluster2 = clusters[j];
0181
0182
0183 auto dz = (static_cast<int>(cluster2.depth()) - static_cast<int>(cluster1.depth()));
0184
0185
0186 if (dz <= 0)
0187 continue;
0188
0189 auto const& crep1 = cluster1.positionREP();
0190 auto const& crep2 = cluster2.positionREP();
0191
0192 auto deta = crep1.eta() - crep2.eta();
0193 deta = deta * deta / (etaRMS2[i] + etaRMS2[j]);
0194 auto dphi = deltaPhi(crep1.phi(), crep2.phi());
0195 dphi = dphi * dphi / (phiRMS2[i] + phiRMS2[j]);
0196
0197
0198
0199 if ((deta < nSigmaEta_) & (dphi < nSigmaPhi_))
0200 links.push_back(ClusterLink(i, j, deta + dphi, std::abs(dz), cluster1.energy() + cluster2.energy()));
0201 }
0202
0203 return links;
0204 }
0205
0206 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::prune(std::vector<ClusterLink>& links,
0207 std::vector<bool>& linkedClusters) {
0208 std::vector<ClusterLink> goodLinks;
0209 std::vector<bool> mask(links.size(), false);
0210 if (links.empty())
0211 return goodLinks;
0212
0213 for (unsigned int i = 0; i < links.size() - 1; ++i) {
0214 if (mask[i])
0215 continue;
0216 for (unsigned int j = i + 1; j < links.size(); ++j) {
0217 if (mask[j])
0218 continue;
0219
0220 const ClusterLink& link1 = links[i];
0221 const ClusterLink& link2 = links[j];
0222
0223 if (link1.to() == link2.to()) {
0224
0225 if (link1.dZ() < link2.dZ()) {
0226 mask[j] = true;
0227 } else if (link1.dZ() > link2.dZ()) {
0228 mask[i] = true;
0229 } else {
0230 if (link1.dR() < link2.dR()) {
0231 mask[j] = true;
0232 } else if (link1.dR() > link2.dR()) {
0233 mask[i] = true;
0234 } else {
0235
0236 if (link1.energy() < link2.energy())
0237 mask[i] = true;
0238 else
0239 mask[j] = true;
0240 }
0241 }
0242 }
0243 }
0244 }
0245
0246 for (unsigned int i = 0; i < links.size(); ++i) {
0247 if (mask[i])
0248 continue;
0249 goodLinks.push_back(links[i]);
0250 linkedClusters[links[i].from()] = true;
0251 linkedClusters[links[i].to()] = true;
0252 }
0253
0254 return goodLinks;
0255 }
0256
0257 void PFMultiDepthClusterizer::absorbCluster(reco::PFCluster& main, const reco::PFCluster& added) {
0258 double e1 = 0.0;
0259 double e2 = 0.0;
0260
0261
0262 for (const auto& fraction : main.recHitFractions())
0263 if (fraction.recHitRef()->detId() == main.seed()) {
0264 e1 = fraction.recHitRef()->energy();
0265 }
0266
0267 for (const auto& fraction : added.recHitFractions()) {
0268 main.addRecHitFraction(fraction);
0269 if (fraction.recHitRef()->detId() == added.seed()) {
0270 e2 = fraction.recHitRef()->energy();
0271 }
0272 }
0273 if (e2 > e1)
0274 main.setSeed(added.seed());
0275 }
0276
0277 void PFMultiDepthClusterizer::expandCluster(reco::PFCluster& cluster,
0278 unsigned int point,
0279 std::vector<bool>& mask,
0280 const reco::PFClusterCollection& clusters,
0281 const std::vector<ClusterLink>& links) {
0282 for (const auto& link : links) {
0283 if (link.from() == point) {
0284
0285 if (!mask[link.from()]) {
0286 absorbCluster(cluster, clusters[link.from()]);
0287 mask[link.from()] = true;
0288 }
0289
0290 if (!mask[link.to()]) {
0291 absorbCluster(cluster, clusters[link.to()]);
0292 mask[link.to()] = true;
0293 expandCluster(cluster, link.to(), mask, clusters, links);
0294 }
0295 }
0296 if (link.to() == point) {
0297
0298 if (!mask[link.to()]) {
0299 absorbCluster(cluster, clusters[link.to()]);
0300 mask[link.to()] = true;
0301 }
0302
0303 if (!mask[link.from()]) {
0304 absorbCluster(cluster, clusters[link.from()]);
0305 mask[link.from()] = true;
0306 expandCluster(cluster, link.from(), mask, clusters, links);
0307 }
0308 }
0309 }
0310 }