PFClusterProducer/plugins/Basic2DGenericPFlowClusterizer.cc

0001 #include "DataFormats/Math/interface/deltaR.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterBuilderBase.h"
0007
0008 #include "Math/GenVector/VectorUtil.h"
0009 #include "vdt/vdtMath.h"
0010
0011 #include <iterator>
0012 #include <unordered_map>
0013
0014 class Basic2DGenericPFlowClusterizer : public PFClusterBuilderBase {
0015   typedef Basic2DGenericPFlowClusterizer B2DGPF;
0016
0017 public:
0018   Basic2DGenericPFlowClusterizer(const edm::ParameterSet& conf, edm::ConsumesCollector& cc);
0019
0020   ~Basic2DGenericPFlowClusterizer() override = default;
0021   Basic2DGenericPFlowClusterizer(const B2DGPF&) = delete;
0022   B2DGPF& operator=(const B2DGPF&) = delete;
0023
0024   void update(const edm::EventSetup& es) override {
0025     _positionCalc->update(es);
0026     if (_allCellsPosCalc)
0027       _allCellsPosCalc->update(es);
0028     if (_convergencePosCalc)
0029       _convergencePosCalc->update(es);
0030   }
0031
0032   void buildClusters(const reco::PFClusterCollection&,
0033                      const std::vector<bool>&,
0034                      reco::PFClusterCollection& outclus) override;
0035
0036 private:
0037   const unsigned _maxIterations;
0038   const double _stoppingTolerance;
0039   const double _showerSigma2;
0040   const bool _excludeOtherSeeds;
0041   const double _minFracTot;
0042   const std::unordered_map<std::string, int> _layerMap;
0043
0044   std::unordered_map<int, std::pair<std::vector<int>, std::vector<double> > > _recHitEnergyNorms;
0045   std::unique_ptr<PFCPositionCalculatorBase> _allCellsPosCalc;
0046   std::unique_ptr<PFCPositionCalculatorBase> _convergencePosCalc;
0047
0048   void seedPFClustersFromTopo(const reco::PFCluster&, const std::vector<bool>&, reco::PFClusterCollection&) const;
0049
0050   void growPFClusters(const reco::PFCluster&,
0051                       const std::vector<bool>&,
0052                       const unsigned toleranceScaling,
0053                       const unsigned iter,
0054                       double dist,
0055                       reco::PFClusterCollection&) const;
0056
0057   void prunePFClusters(reco::PFClusterCollection&) const;
0058 };
0059
0060 DEFINE_EDM_PLUGIN(PFClusterBuilderFactory, Basic2DGenericPFlowClusterizer, "Basic2DGenericPFlowClusterizer");
0061
0062 #ifdef PFLOW_DEBUG
0063 #define LOGVERB(x) edm::LogVerbatim(x)
0064 #define LOGWARN(x) edm::LogWarning(x)
0065 #define LOGERR(x) edm::LogError(x)
0066 #define LOGDRESSED(x) edm::LogInfo(x)
0067 #else
0068 #define LOGVERB(x) LogTrace(x)
0069 #define LOGWARN(x) edm::LogWarning(x)
0070 #define LOGERR(x) edm::LogError(x)
0071 #define LOGDRESSED(x) LogDebug(x)
0072 #endif
0073
0074 Basic2DGenericPFlowClusterizer::Basic2DGenericPFlowClusterizer(const edm::ParameterSet& conf,
0075                                                                edm::ConsumesCollector& cc)
0076     : PFClusterBuilderBase(conf, cc),
0077       _maxIterations(conf.getParameter<unsigned>("maxIterations")),
0078       _stoppingTolerance(conf.getParameter<double>("stoppingTolerance")),
0079       _showerSigma2(std::pow(conf.getParameter<double>("showerSigma"), 2.0)),
0080       _excludeOtherSeeds(conf.getParameter<bool>("excludeOtherSeeds")),
0081       _minFracTot(conf.getParameter<double>("minFracTot")),
0082       _layerMap({{"PS2", (int)PFLayer::PS2},
0083                  {"PS1", (int)PFLayer::PS1},
0084                  {"ECAL_ENDCAP", (int)PFLayer::ECAL_ENDCAP},
0085                  {"ECAL_BARREL", (int)PFLayer::ECAL_BARREL},
0086                  {"NONE", (int)PFLayer::NONE},
0087                  {"HCAL_BARREL1", (int)PFLayer::HCAL_BARREL1},
0088                  {"HCAL_BARREL2_RING0", (int)PFLayer::HCAL_BARREL2},
0089                  {"HCAL_BARREL2_RING1", 100 * (int)PFLayer::HCAL_BARREL2},
0090                  {"HCAL_ENDCAP", (int)PFLayer::HCAL_ENDCAP},
0091                  {"HF_EM", (int)PFLayer::HF_EM},
0092                  {"HF_HAD", (int)PFLayer::HF_HAD}}) {
0093   const std::vector<edm::ParameterSet>& thresholds = conf.getParameterSetVector("recHitEnergyNorms");
0094   for (const auto& pset : thresholds) {
0095     const std::string& det = pset.getParameter<std::string>("detector");
0096
0097     std::vector<int> depths;
0098     std::vector<double> rhE_norm;
0099
0100     if (det == std::string("HCAL_BARREL1") || det == std::string("HCAL_ENDCAP")) {
0101       depths = pset.getParameter<std::vector<int> >("depths");
0102       rhE_norm = pset.getParameter<std::vector<double> >("recHitEnergyNorm");
0103     } else {
0104       depths.push_back(0);
0105       rhE_norm.push_back(pset.getParameter<double>("recHitEnergyNorm"));
0106     }
0107
0108     if (rhE_norm.size() != depths.size()) {
0109       throw cms::Exception("InvalidPFRecHitThreshold")
0110           << "PFlowClusterizerThreshold mismatch with the numbers of depths";
0111     }
0112
0113     auto entry = _layerMap.find(det);
0114     if (entry == _layerMap.end()) {
0115       throw cms::Exception("InvalidDetectorLayer") << "Detector layer : " << det << " is not in the list of recognized"
0116                                                    << " detector layers!";
0117     }
0118     _recHitEnergyNorms.emplace(_layerMap.find(det)->second, std::make_pair(depths, rhE_norm));
0119   }
0120
0121   if (conf.exists("allCellsPositionCalc")) {
0122     const edm::ParameterSet& acConf = conf.getParameterSet("allCellsPositionCalc");
0123     const std::string& algoac = acConf.getParameter<std::string>("algoName");
0124     _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0125   }
0126   // if necessary a third pos calc for convergence testing
0127   if (conf.exists("positionCalcForConvergence")) {
0128     const edm::ParameterSet& convConf = conf.getParameterSet("positionCalcForConvergence");
0129     const std::string& algoconv = convConf.getParameter<std::string>("algoName");
0130     _convergencePosCalc = PFCPositionCalculatorFactory::get()->create(algoconv, convConf, cc);
0131   }
0132 }
0133
0134 void Basic2DGenericPFlowClusterizer::buildClusters(const reco::PFClusterCollection& input,
0135                                                    const std::vector<bool>& seedable,
0136                                                    reco::PFClusterCollection& output) {
0137   reco::PFClusterCollection clustersInTopo;
0138   for (const auto& topocluster : input) {
0139     clustersInTopo.clear();
0140     seedPFClustersFromTopo(topocluster, seedable, clustersInTopo);
0141     const unsigned tolScal = std::pow(std::max(1.0, clustersInTopo.size() - 1.0), 2.0);
0142     growPFClusters(topocluster, seedable, tolScal, 0, tolScal, clustersInTopo);
0143     // step added by Josh Bendavid, removes low-fraction clusters
0144     // did not impact position resolution with fraction cut of 1e-7
0145     // decreases the size of each pf cluster considerably
0146     prunePFClusters(clustersInTopo);
0147     // recalculate the positions of the pruned clusters
0148     if (_convergencePosCalc) {
0149       // if defined, use the special position calculation for convergence tests
0150       _convergencePosCalc->calculateAndSetPositions(clustersInTopo);
0151     } else {
0152       if (clustersInTopo.size() == 1 && _allCellsPosCalc) {
0153         _allCellsPosCalc->calculateAndSetPosition(clustersInTopo.back());
0154       } else {
0155         _positionCalc->calculateAndSetPositions(clustersInTopo);
0156       }
0157     }
0158     for (auto& clusterout : clustersInTopo) {
0159       output.insert(output.end(), std::move(clusterout));
0160     }
0161   }
0162 }
0163
0164 void Basic2DGenericPFlowClusterizer::seedPFClustersFromTopo(const reco::PFCluster& topo,
0165                                                             const std::vector<bool>& seedable,
0166                                                             reco::PFClusterCollection& initialPFClusters) const {
0167   const auto& recHitFractions = topo.recHitFractions();
0168   for (const auto& rhf : recHitFractions) {
0169     if (!seedable[rhf.recHitRef().key()])
0170       continue;
0171     initialPFClusters.push_back(reco::PFCluster());
0172     reco::PFCluster& current = initialPFClusters.back();
0173     current.addRecHitFraction(rhf);
0174     current.setSeed(rhf.recHitRef()->detId());
0175     if (_convergencePosCalc) {
0176       _convergencePosCalc->calculateAndSetPosition(current);
0177     } else {
0178       _positionCalc->calculateAndSetPosition(current);
0179     }
0180   }
0181 }
0182
0183 void Basic2DGenericPFlowClusterizer::growPFClusters(const reco::PFCluster& topo,
0184                                                     const std::vector<bool>& seedable,
0185                                                     const unsigned toleranceScaling,
0186                                                     const unsigned iter,
0187                                                     double diff,
0188                                                     reco::PFClusterCollection& clusters) const {
0189   if (iter >= _maxIterations) {
0190     LOGDRESSED("Basic2DGenericPFlowClusterizer:growAndStabilizePFClusters")
0191         << "reached " << _maxIterations << " iterations, terminated position "
0192         << "fit with diff = " << diff;
0193   }
0194   if (iter >= _maxIterations || diff <= _stoppingTolerance * toleranceScaling)
0195     return;
0196   // reset the rechits in this cluster, keeping the previous position
0197   std::vector<reco::PFCluster::REPPoint> clus_prev_pos;
0198   for (auto& cluster : clusters) {
0199     const reco::PFCluster::REPPoint& repp = cluster.positionREP();
0200     clus_prev_pos.emplace_back(repp.rho(), repp.eta(), repp.phi());
0201     if (_convergencePosCalc) {
0202       if (clusters.size() == 1 && _allCellsPosCalc) {
0203         _allCellsPosCalc->calculateAndSetPosition(cluster);
0204       } else {
0205         _positionCalc->calculateAndSetPosition(cluster);
0206       }
0207     }
0208     cluster.resetHitsAndFractions();
0209   }
0210   // loop over topo cluster and grow current PFCluster hypothesis
0211   std::vector<double> dist2, frac;
0212   double fractot = 0;
0213   for (const reco::PFRecHitFraction& rhf : topo.recHitFractions()) {
0214     const reco::PFRecHitRef& refhit = rhf.recHitRef();
0215     int cell_layer = (int)refhit->layer();
0216     if (cell_layer == PFLayer::HCAL_BARREL2 && std::abs(refhit->positionREP().eta()) > 0.34) {
0217       cell_layer *= 100;
0218     }
0219
0220     math::XYZPoint topocellpos_xyz(refhit->position());
0221     dist2.clear();
0222     frac.clear();
0223     fractot = 0;
0224
0225     double recHitEnergyNorm = 0.;
0226     auto const& recHitEnergyNormDepthPair = _recHitEnergyNorms.find(cell_layer)->second;
0227
0228     for (unsigned int j = 0; j < recHitEnergyNormDepthPair.second.size(); ++j) {
0229       int depth = recHitEnergyNormDepthPair.first[j];
0230
0231       if ((cell_layer == PFLayer::HCAL_BARREL1 && refhit->depth() == depth) ||
0232           (cell_layer == PFLayer::HCAL_ENDCAP && refhit->depth() == depth) ||
0233           (cell_layer != PFLayer::HCAL_ENDCAP && cell_layer != PFLayer::HCAL_BARREL1))
0234         recHitEnergyNorm = recHitEnergyNormDepthPair.second[j];
0235     }
0236
0237     // add rechits to clusters, calculating fraction based on distance
0238     for (auto& cluster : clusters) {
0239       const math::XYZPoint& clusterpos_xyz = cluster.position();
0240       const math::XYZVector deltav = clusterpos_xyz - topocellpos_xyz;
0241       const double d2 = deltav.Mag2() / _showerSigma2;
0242       dist2.emplace_back(d2);
0243       if (d2 > 100) {
0244         LOGDRESSED("Basic2DGenericPFlowClusterizer:growAndStabilizePFClusters")
0245             << "Warning! :: pfcluster-topocell distance is too large! d= " << d2;
0246       }
0247
0248       // fraction assignment logic
0249       double fraction;
0250       if (refhit->detId() == cluster.seed() && _excludeOtherSeeds) {
0251         fraction = 1.0;
0252       } else if (seedable[refhit.key()] && _excludeOtherSeeds) {
0253         fraction = 0.0;
0254       } else {
0255         fraction = cluster.energy() / recHitEnergyNorm * vdt::fast_expf(-0.5 * d2);
0256       }
0257       fractot += fraction;
0258       frac.emplace_back(fraction);
0259     }
0260     for (unsigned i = 0; i < clusters.size(); ++i) {
0261       if (fractot > _minFracTot || (refhit->detId() == clusters[i].seed() && fractot > 0.0)) {
0262         frac[i] /= fractot;
0263       } else {
0264         continue;
0265       }
0266       // if the fraction has been set to 0, the cell
0267       // is now added to the cluster - careful ! (PJ, 19/07/08)
0268       // BUT KEEP ONLY CLOSE CELLS OTHERWISE MEMORY JUST EXPLOSES
0269       // (PJ, 15/09/08 <- similar to what existed before the
0270       // previous bug fix, but keeps the close seeds inside,
0271       // even if their fraction was set to zero.)
0272       // Also add a protection to keep the seed in the cluster
0273       // when the latter gets far from the former. These cases
0274       // (about 1% of the clusters) need to be studied, as
0275       // they create fake photons, in general.
0276       // (PJ, 16/09/08)
0277       if (dist2[i] < 100.0 || frac[i] > 0.9999) {
0278         clusters[i].addRecHitFraction(reco::PFRecHitFraction(refhit, frac[i]));
0279       }
0280     }
0281   }
0282   // recalculate positions and calculate convergence parameter
0283   double diff2 = 0.0;
0284   for (unsigned i = 0; i < clusters.size(); ++i) {
0285     if (_convergencePosCalc) {
0286       _convergencePosCalc->calculateAndSetPosition(clusters[i]);
0287     } else {
0288       if (clusters.size() == 1 && _allCellsPosCalc) {
0289         _allCellsPosCalc->calculateAndSetPosition(clusters[i]);
0290       } else {
0291         _positionCalc->calculateAndSetPosition(clusters[i]);
0292       }
0293     }
0294     const double delta2 = reco::deltaR2(clusters[i].positionREP(), clus_prev_pos[i]);
0295     if (delta2 > diff2)
0296       diff2 = delta2;
0297   }
0298   diff = std::sqrt(diff2);
0299   dist2.clear();
0300   frac.clear();
0301   clus_prev_pos.clear();  // avoid badness
0302   growPFClusters(topo, seedable, toleranceScaling, iter + 1, diff, clusters);
0303 }
0304
0305 void Basic2DGenericPFlowClusterizer::prunePFClusters(reco::PFClusterCollection& clusters) const {
0306   for (auto& cluster : clusters) {
0307     cluster.pruneUsing([&](const reco::PFRecHitFraction& rhf) { return rhf.fraction() > _minFractionToKeep; });
0308   }
0309 }