Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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