PFClusterProducer/plugins/PFMultiDepthClusterizer.cc

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
0010 #include "vdt/vdtMath.h"
0011
0012 #include <iterator>
0013 #include <unordered_map>
0014
0015 class PFMultiDepthClusterizer final : public PFClusterBuilderBase {
0016   typedef PFMultiDepthClusterizer B2DGPF;
0017
0018 public:
0019   PFMultiDepthClusterizer(const edm::ParameterSet& conf, edm::ConsumesCollector& cc);
0020
0021   ~PFMultiDepthClusterizer() override = default;
0022   PFMultiDepthClusterizer(const B2DGPF&) = delete;
0023   B2DGPF& operator=(const B2DGPF&) = delete;
0024
0025   void update(const edm::EventSetup& es) override { _allCellsPosCalc->update(es); }
0026
0027   void buildClusters(const reco::PFClusterCollection&,
0028                      const std::vector<bool>&,
0029                      reco::PFClusterCollection& outclus) override;
0030
0031 private:
0032   std::unique_ptr<PFCPositionCalculatorBase> _allCellsPosCalc;
0033   double nSigmaEta_;
0034   double nSigmaPhi_;
0035
0036   class ClusterLink {
0037   public:
0038     ClusterLink(unsigned int i, unsigned int j, double DR, int DZ, double energy) {
0039       from_ = i;
0040       to_ = j;
0041       linkDR_ = DR;
0042       linkDZ_ = DZ;
0043       linkE_ = energy;
0044     }
0045
0046     ~ClusterLink() = default;
0047
0048     unsigned int from() const { return from_; }
0049     unsigned int to() const { return to_; }
0050     double dR() const { return linkDR_; }
0051     int dZ() const { return linkDZ_; }
0052     double energy() const { return linkE_; }
0053
0054   private:
0055     unsigned int from_;
0056     unsigned int to_;
0057     double linkDR_;
0058     int linkDZ_;
0059     double linkE_;
0060   };
0061
0062   void calculateShowerShapes(const reco::PFClusterCollection&, std::vector<double>&, std::vector<double>&);
0063   std::vector<ClusterLink> link(const reco::PFClusterCollection&,
0064                                 const std::vector<double>&,
0065                                 const std::vector<double>&);
0066   std::vector<ClusterLink> prune(std::vector<ClusterLink>&, std::vector<bool>& linkedClusters);
0067
0068   void expandCluster(reco::PFCluster&,
0069                      unsigned int point,
0070                      std::vector<bool>& mask,
0071                      const reco::PFClusterCollection&,
0072                      const std::vector<ClusterLink>& links);
0073
0074   void absorbCluster(reco::PFCluster&, const reco::PFCluster&);
0075 };
0076
0077 DEFINE_EDM_PLUGIN(PFClusterBuilderFactory, PFMultiDepthClusterizer, "PFMultiDepthClusterizer");
0078
0079 PFMultiDepthClusterizer::PFMultiDepthClusterizer(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0080     : PFClusterBuilderBase(conf, cc) {
0081   if (conf.exists("allCellsPositionCalc")) {
0082     const edm::ParameterSet& acConf = conf.getParameterSet("allCellsPositionCalc");
0083     const std::string& algoac = acConf.getParameter<std::string>("algoName");
0084     _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0085   }
0086
0087   nSigmaEta_ = pow(conf.getParameter<double>("nSigmaEta"), 2);
0088   nSigmaPhi_ = pow(conf.getParameter<double>("nSigmaPhi"), 2);
0089 }
0090
0091 void PFMultiDepthClusterizer::buildClusters(const reco::PFClusterCollection& input,
0092                                             const std::vector<bool>& seedable,
0093                                             reco::PFClusterCollection& output) {
0094   std::vector<double> etaRMS2(input.size(), 0.0);
0095   std::vector<double> phiRMS2(input.size(), 0.0);
0096
0097   //We need to sort the clusters for smaller to larger depth
0098   //  for (unsigned int i=0;i<input.size();++i)
0099   //   printf(" cluster%f %f \n",input[i].depth(),input[i].energy());
0100
0101   //calculate cluster shapes
0102   calculateShowerShapes(input, etaRMS2, phiRMS2);
0103
0104   //link
0105   auto&& links = link(input, etaRMS2, phiRMS2);
0106   //  for (const auto& link: links)
0107   //    printf("link %d %d %f %f\n",link.from(),link.to(),link.dR(),link.dZ());
0108
0109   std::vector<bool> mask(input.size(), false);
0110   std::vector<bool> linked(input.size(), false);
0111
0112   //prune
0113   auto&& prunedLinks = prune(links, linked);
0114
0115   //printf("Pruned links\n")
0116   //  for (const auto& link: prunedLinks)
0117   //  printf("link %d %d %f %f\n",link.from(),link.to(),link.dR(),link.dZ());
0118
0119   //now we need to build clusters
0120   for (unsigned int i = 0; i < input.size(); ++i) {
0121     //if not masked
0122     if (mask[i])
0123       continue;
0124     //if not linked just spit it out
0125     if (!linked[i]) {
0126       output.push_back(input[i]);
0127       //      printf("Added single cluster with energy =%f \n",input[i].energy());
0128       mask[i] = true;
0129       continue;
0130     }
0131
0132     //now business: if  linked and not  masked gather clusters
0133     reco::PFCluster cluster = input[i];
0134     mask[i] = true;
0135     expandCluster(cluster, i, mask, input, prunedLinks);
0136     _allCellsPosCalc->calculateAndSetPosition(cluster);
0137     output.push_back(cluster);
0138     //    printf("Added linked cluster with energy =%f\n",cluster.energy());
0139   }
0140 }
0141
0142 void PFMultiDepthClusterizer::calculateShowerShapes(const reco::PFClusterCollection& clusters,
0143                                                     std::vector<double>& etaRMS2,
0144                                                     std::vector<double>& phiRMS2) {
0145   //shower shapes. here do not use the fractions
0146
0147   for (unsigned int i = 0; i < clusters.size(); ++i) {
0148     const reco::PFCluster& cluster = clusters[i];
0149     double etaSum = 0.0;
0150     double phiSum = 0.0;
0151     auto const& crep = cluster.positionREP();
0152     for (const auto& frac : cluster.recHitFractions()) {
0153       auto const& h = *frac.recHitRef();
0154       auto const& rep = h.positionREP();
0155       etaSum += (frac.fraction() * h.energy()) * std::abs(rep.eta() - crep.eta());
0156       phiSum += (frac.fraction() * h.energy()) * std::abs(deltaPhi(rep.phi(), crep.phi()));
0157     }
0158     //protection for single line : assign ~ tower
0159     etaRMS2[i] = std::max(etaSum / cluster.energy(), 0.1);
0160     etaRMS2[i] *= etaRMS2[i];
0161     phiRMS2[i] = std::max(phiSum / cluster.energy(), 0.1);
0162     phiRMS2[i] *= phiRMS2[i];
0163   }
0164 }
0165
0166 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::link(
0167     const reco::PFClusterCollection& clusters, const std::vector<double>& etaRMS2, const std::vector<double>& phiRMS2) {
0168   std::vector<ClusterLink> links;
0169   //loop on all pairs
0170   for (unsigned int i = 0; i < clusters.size(); ++i)
0171     for (unsigned int j = 0; j < clusters.size(); ++j) {
0172       if (i == j)
0173         continue;
0174
0175       const reco::PFCluster& cluster1 = clusters[i];
0176       const reco::PFCluster& cluster2 = clusters[j];
0177
0178       // PFCluster depth stored as double but HCAL layer clusters have integral depths only
0179       auto dz = (static_cast<int>(cluster2.depth()) - static_cast<int>(cluster1.depth()));
0180
0181       //Do not link at the same layer and only link inside out!
0182       if (dz <= 0)
0183         continue;
0184
0185       auto const& crep1 = cluster1.positionREP();
0186       auto const& crep2 = cluster2.positionREP();
0187
0188       auto deta = crep1.eta() - crep2.eta();
0189       deta = deta * deta / (etaRMS2[i] + etaRMS2[j]);
0190       auto dphi = deltaPhi(crep1.phi(), crep2.phi());
0191       dphi = dphi * dphi / (phiRMS2[i] + phiRMS2[j]);
0192
0193       //      printf("Testing Link %d -> %d (%f %f %f %f ) \n",i,j,deta,dphi,cluster1.position().Eta()-cluster2.position().Eta(),deltaPhi(cluster1.position().Phi(),cluster2.position().Phi()));
0194
0195       if ((deta < nSigmaEta_) & (dphi < nSigmaPhi_))
0196         links.push_back(ClusterLink(i, j, deta + dphi, std::abs(dz), cluster1.energy() + cluster2.energy()));
0197     }
0198
0199   return links;
0200 }
0201
0202 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::prune(std::vector<ClusterLink>& links,
0203                                                                                  std::vector<bool>& linkedClusters) {
0204   std::vector<ClusterLink> goodLinks;
0205   std::vector<bool> mask(links.size(), false);
0206   if (links.empty())
0207     return goodLinks;
0208
0209   for (unsigned int i = 0; i < links.size() - 1; ++i) {
0210     if (mask[i])
0211       continue;
0212     for (unsigned int j = i + 1; j < links.size(); ++j) {
0213       if (mask[j])
0214         continue;
0215
0216       const ClusterLink& link1 = links[i];
0217       const ClusterLink& link2 = links[j];
0218
0219       if (link1.to() == link2.to()) {  //found two links going to the same spot,kill one
0220         //first prefer nearby layers
0221         if (link1.dZ() < link2.dZ()) {
0222           mask[j] = true;
0223         } else if (link1.dZ() > link2.dZ()) {
0224           mask[i] = true;
0225         } else {  //if same layer-pick based on transverse compatibility
0226           if (link1.dR() < link2.dR()) {
0227             mask[j] = true;
0228           } else if (link1.dR() > link2.dR()) {
0229             mask[i] = true;
0230           } else {
0231             //same distance as well -> can happen !!!!! Pick the highest SUME
0232             if (link1.energy() < link2.energy())
0233               mask[i] = true;
0234             else
0235               mask[j] = true;
0236           }
0237         }
0238       }
0239     }
0240   }
0241
0242   for (unsigned int i = 0; i < links.size(); ++i) {
0243     if (mask[i])
0244       continue;
0245     goodLinks.push_back(links[i]);
0246     linkedClusters[links[i].from()] = true;
0247     linkedClusters[links[i].to()] = true;
0248   }
0249
0250   return goodLinks;
0251 }
0252
0253 void PFMultiDepthClusterizer::absorbCluster(reco::PFCluster& main, const reco::PFCluster& added) {
0254   double e1 = 0.0;
0255   double e2 = 0.0;
0256
0257   //find seeds
0258   for (const auto& fraction : main.recHitFractions())
0259     if (fraction.recHitRef()->detId() == main.seed()) {
0260       e1 = fraction.recHitRef()->energy();
0261     }
0262
0263   for (const auto& fraction : added.recHitFractions()) {
0264     main.addRecHitFraction(fraction);
0265     if (fraction.recHitRef()->detId() == added.seed()) {
0266       e2 = fraction.recHitRef()->energy();
0267     }
0268   }
0269   if (e2 > e1)
0270     main.setSeed(added.seed());
0271 }
0272
0273 void PFMultiDepthClusterizer::expandCluster(reco::PFCluster& cluster,
0274                                             unsigned int point,
0275                                             std::vector<bool>& mask,
0276                                             const reco::PFClusterCollection& clusters,
0277                                             const std::vector<ClusterLink>& links) {
0278   for (const auto& link : links) {
0279     if (link.from() == point) {
0280       //found link that starts from this guy if not masked absorb
0281       if (!mask[link.from()]) {
0282         absorbCluster(cluster, clusters[link.from()]);
0283         mask[link.from()] = true;
0284       }
0285
0286       if (!mask[link.to()]) {
0287         absorbCluster(cluster, clusters[link.to()]);
0288         mask[link.to()] = true;
0289         expandCluster(cluster, link.to(), mask, clusters, links);
0290       }
0291     }
0292     if (link.to() == point) {
0293       //found link that starts from this guy if not masked absorb
0294       if (!mask[link.to()]) {
0295         absorbCluster(cluster, clusters[link.to()]);
0296         mask[link.to()] = true;
0297       }
0298
0299       if (!mask[link.from()]) {
0300         absorbCluster(cluster, clusters[link.from()]);
0301         mask[link.from()] = true;
0302         expandCluster(cluster, link.from(), mask, clusters, links);
0303       }
0304     }
0305   }
0306 }