Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:27:25

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   if (conf.exists("allCellsPositionCalc")) {
0084     const edm::ParameterSet& acConf = conf.getParameterSet("allCellsPositionCalc");
0085     const std::string& algoac = acConf.getParameter<std::string>("algoName");
0086     _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0087   }
0088 
0089   nSigmaEta_ = pow(conf.getParameter<double>("nSigmaEta"), 2);
0090   nSigmaPhi_ = pow(conf.getParameter<double>("nSigmaPhi"), 2);
0091 }
0092 
0093 void PFMultiDepthClusterizer::buildClusters(const reco::PFClusterCollection& input,
0094                                             const std::vector<bool>& seedable,
0095                                             reco::PFClusterCollection& output,
0096                                             const HcalPFCuts* hcalCuts) {
0097   std::vector<double> etaRMS2(input.size(), 0.0);
0098   std::vector<double> phiRMS2(input.size(), 0.0);
0099 
0100   //We need to sort the clusters for smaller to larger depth
0101   //  for (unsigned int i=0;i<input.size();++i)
0102   //   printf(" cluster%f %f \n",input[i].depth(),input[i].energy());
0103 
0104   //calculate cluster shapes
0105   calculateShowerShapes(input, etaRMS2, phiRMS2);
0106 
0107   //link
0108   auto&& links = link(input, etaRMS2, phiRMS2);
0109   //  for (const auto& link: links)
0110   //    printf("link %d %d %f %f\n",link.from(),link.to(),link.dR(),link.dZ());
0111 
0112   std::vector<bool> mask(input.size(), false);
0113   std::vector<bool> linked(input.size(), false);
0114 
0115   //prune
0116   auto&& prunedLinks = prune(links, linked);
0117 
0118   //printf("Pruned links\n")
0119   //  for (const auto& link: prunedLinks)
0120   //  printf("link %d %d %f %f\n",link.from(),link.to(),link.dR(),link.dZ());
0121 
0122   //now we need to build clusters
0123   for (unsigned int i = 0; i < input.size(); ++i) {
0124     //if not masked
0125     if (mask[i])
0126       continue;
0127     //if not linked just spit it out
0128     if (!linked[i]) {
0129       output.push_back(input[i]);
0130       //      printf("Added single cluster with energy =%f \n",input[i].energy());
0131       mask[i] = true;
0132       continue;
0133     }
0134 
0135     //now business: if  linked and not  masked gather clusters
0136     reco::PFCluster cluster = input[i];
0137     mask[i] = true;
0138     expandCluster(cluster, i, mask, input, prunedLinks);
0139     _allCellsPosCalc->calculateAndSetPosition(cluster, hcalCuts);
0140     output.push_back(cluster);
0141     //    printf("Added linked cluster with energy =%f\n",cluster.energy());
0142   }
0143 }
0144 
0145 void PFMultiDepthClusterizer::calculateShowerShapes(const reco::PFClusterCollection& clusters,
0146                                                     std::vector<double>& etaRMS2,
0147                                                     std::vector<double>& phiRMS2) {
0148   //shower shapes. here do not use the fractions
0149 
0150   for (unsigned int i = 0; i < clusters.size(); ++i) {
0151     const reco::PFCluster& cluster = clusters[i];
0152     double etaSum = 0.0;
0153     double phiSum = 0.0;
0154     auto const& crep = cluster.positionREP();
0155     for (const auto& frac : cluster.recHitFractions()) {
0156       auto const& h = *frac.recHitRef();
0157       auto const& rep = h.positionREP();
0158       etaSum += (frac.fraction() * h.energy()) * std::abs(rep.eta() - crep.eta());
0159       phiSum += (frac.fraction() * h.energy()) * std::abs(deltaPhi(rep.phi(), crep.phi()));
0160     }
0161     //protection for single line : assign ~ tower
0162     etaRMS2[i] = std::max(etaSum / cluster.energy(), 0.1);
0163     etaRMS2[i] *= etaRMS2[i];
0164     phiRMS2[i] = std::max(phiSum / cluster.energy(), 0.1);
0165     phiRMS2[i] *= phiRMS2[i];
0166   }
0167 }
0168 
0169 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::link(
0170     const reco::PFClusterCollection& clusters, const std::vector<double>& etaRMS2, const std::vector<double>& phiRMS2) {
0171   std::vector<ClusterLink> links;
0172   //loop on all pairs
0173   for (unsigned int i = 0; i < clusters.size(); ++i)
0174     for (unsigned int j = 0; j < clusters.size(); ++j) {
0175       if (i == j)
0176         continue;
0177 
0178       const reco::PFCluster& cluster1 = clusters[i];
0179       const reco::PFCluster& cluster2 = clusters[j];
0180 
0181       // PFCluster depth stored as double but HCAL layer clusters have integral depths only
0182       auto dz = (static_cast<int>(cluster2.depth()) - static_cast<int>(cluster1.depth()));
0183 
0184       //Do not link at the same layer and only link inside out!
0185       if (dz <= 0)
0186         continue;
0187 
0188       auto const& crep1 = cluster1.positionREP();
0189       auto const& crep2 = cluster2.positionREP();
0190 
0191       auto deta = crep1.eta() - crep2.eta();
0192       deta = deta * deta / (etaRMS2[i] + etaRMS2[j]);
0193       auto dphi = deltaPhi(crep1.phi(), crep2.phi());
0194       dphi = dphi * dphi / (phiRMS2[i] + phiRMS2[j]);
0195 
0196       //      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()));
0197 
0198       if ((deta < nSigmaEta_) & (dphi < nSigmaPhi_))
0199         links.push_back(ClusterLink(i, j, deta + dphi, std::abs(dz), cluster1.energy() + cluster2.energy()));
0200     }
0201 
0202   return links;
0203 }
0204 
0205 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::prune(std::vector<ClusterLink>& links,
0206                                                                                  std::vector<bool>& linkedClusters) {
0207   std::vector<ClusterLink> goodLinks;
0208   std::vector<bool> mask(links.size(), false);
0209   if (links.empty())
0210     return goodLinks;
0211 
0212   for (unsigned int i = 0; i < links.size() - 1; ++i) {
0213     if (mask[i])
0214       continue;
0215     for (unsigned int j = i + 1; j < links.size(); ++j) {
0216       if (mask[j])
0217         continue;
0218 
0219       const ClusterLink& link1 = links[i];
0220       const ClusterLink& link2 = links[j];
0221 
0222       if (link1.to() == link2.to()) {  //found two links going to the same spot,kill one
0223         //first prefer nearby layers
0224         if (link1.dZ() < link2.dZ()) {
0225           mask[j] = true;
0226         } else if (link1.dZ() > link2.dZ()) {
0227           mask[i] = true;
0228         } else {  //if same layer-pick based on transverse compatibility
0229           if (link1.dR() < link2.dR()) {
0230             mask[j] = true;
0231           } else if (link1.dR() > link2.dR()) {
0232             mask[i] = true;
0233           } else {
0234             //same distance as well -> can happen !!!!! Pick the highest SUME
0235             if (link1.energy() < link2.energy())
0236               mask[i] = true;
0237             else
0238               mask[j] = true;
0239           }
0240         }
0241       }
0242     }
0243   }
0244 
0245   for (unsigned int i = 0; i < links.size(); ++i) {
0246     if (mask[i])
0247       continue;
0248     goodLinks.push_back(links[i]);
0249     linkedClusters[links[i].from()] = true;
0250     linkedClusters[links[i].to()] = true;
0251   }
0252 
0253   return goodLinks;
0254 }
0255 
0256 void PFMultiDepthClusterizer::absorbCluster(reco::PFCluster& main, const reco::PFCluster& added) {
0257   double e1 = 0.0;
0258   double e2 = 0.0;
0259 
0260   //find seeds
0261   for (const auto& fraction : main.recHitFractions())
0262     if (fraction.recHitRef()->detId() == main.seed()) {
0263       e1 = fraction.recHitRef()->energy();
0264     }
0265 
0266   for (const auto& fraction : added.recHitFractions()) {
0267     main.addRecHitFraction(fraction);
0268     if (fraction.recHitRef()->detId() == added.seed()) {
0269       e2 = fraction.recHitRef()->energy();
0270     }
0271   }
0272   if (e2 > e1)
0273     main.setSeed(added.seed());
0274 }
0275 
0276 void PFMultiDepthClusterizer::expandCluster(reco::PFCluster& cluster,
0277                                             unsigned int point,
0278                                             std::vector<bool>& mask,
0279                                             const reco::PFClusterCollection& clusters,
0280                                             const std::vector<ClusterLink>& links) {
0281   for (const auto& link : links) {
0282     if (link.from() == point) {
0283       //found link that starts from this guy if not masked absorb
0284       if (!mask[link.from()]) {
0285         absorbCluster(cluster, clusters[link.from()]);
0286         mask[link.from()] = true;
0287       }
0288 
0289       if (!mask[link.to()]) {
0290         absorbCluster(cluster, clusters[link.to()]);
0291         mask[link.to()] = true;
0292         expandCluster(cluster, link.to(), mask, clusters, links);
0293       }
0294     }
0295     if (link.to() == point) {
0296       //found link that starts from this guy if not masked absorb
0297       if (!mask[link.to()]) {
0298         absorbCluster(cluster, clusters[link.to()]);
0299         mask[link.to()] = true;
0300       }
0301 
0302       if (!mask[link.from()]) {
0303         absorbCluster(cluster, clusters[link.from()]);
0304         mask[link.from()] = true;
0305         expandCluster(cluster, link.from(), mask, clusters, links);
0306       }
0307     }
0308   }
0309 }