File indexing completed on 2024-06-22 02:24:03
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 const auto& acConf = conf.getParameterSet("allCellsPositionCalc");
0129 if (!acConf.empty()) {
0130 const std::string& algoac = acConf.getParameter<std::string>("algoName");
0131 if (!algoac.empty())
0132 _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0133 }
0134
0135 const auto& convConf = conf.getParameterSet("positionCalcForConvergence");
0136 if (!convConf.empty()) {
0137 const std::string& algoconv = convConf.getParameter<std::string>("algoName");
0138 if (!algoconv.empty())
0139 _convergencePosCalc = PFCPositionCalculatorFactory::get()->create(algoconv, convConf, cc);
0140 }
0141 }
0142
0143 void Basic2DGenericPFlowClusterizer::buildClusters(const reco::PFClusterCollection& input,
0144 const std::vector<bool>& seedable,
0145 reco::PFClusterCollection& output,
0146 const HcalPFCuts* hcalCuts) {
0147 reco::PFClusterCollection clustersInTopo;
0148 for (const auto& topocluster : input) {
0149 clustersInTopo.clear();
0150 seedPFClustersFromTopo(topocluster, seedable, clustersInTopo, hcalCuts);
0151 const unsigned tolScal = std::pow(std::max(1.0, clustersInTopo.size() - 1.0), 2.0);
0152 growPFClusters(topocluster, seedable, tolScal, 0, tolScal, clustersInTopo, hcalCuts);
0153
0154
0155
0156 prunePFClusters(clustersInTopo);
0157
0158 if (_convergencePosCalc) {
0159
0160 _convergencePosCalc->calculateAndSetPositions(clustersInTopo, hcalCuts);
0161 } else {
0162 if (clustersInTopo.size() == 1 && _allCellsPosCalc) {
0163 _allCellsPosCalc->calculateAndSetPosition(clustersInTopo.back(), hcalCuts);
0164 } else {
0165 _positionCalc->calculateAndSetPositions(clustersInTopo, hcalCuts);
0166 }
0167 }
0168 for (auto& clusterout : clustersInTopo) {
0169 output.insert(output.end(), std::move(clusterout));
0170 }
0171 }
0172 }
0173
0174 void Basic2DGenericPFlowClusterizer::seedPFClustersFromTopo(const reco::PFCluster& topo,
0175 const std::vector<bool>& seedable,
0176 reco::PFClusterCollection& initialPFClusters,
0177 const HcalPFCuts* hcalCuts) const {
0178 const auto& recHitFractions = topo.recHitFractions();
0179 for (const auto& rhf : recHitFractions) {
0180 if (!seedable[rhf.recHitRef().key()])
0181 continue;
0182 initialPFClusters.push_back(reco::PFCluster());
0183 reco::PFCluster& current = initialPFClusters.back();
0184 current.addRecHitFraction(rhf);
0185 current.setSeed(rhf.recHitRef()->detId());
0186 if (_convergencePosCalc) {
0187 _convergencePosCalc->calculateAndSetPosition(current, hcalCuts);
0188 } else {
0189 _positionCalc->calculateAndSetPosition(current, hcalCuts);
0190 }
0191 }
0192 }
0193
0194 void Basic2DGenericPFlowClusterizer::growPFClusters(const reco::PFCluster& topo,
0195 const std::vector<bool>& seedable,
0196 const unsigned toleranceScaling,
0197 const unsigned iter,
0198 double diff,
0199 reco::PFClusterCollection& clusters,
0200 const HcalPFCuts* hcalCuts) const {
0201 if (iter >= _maxIterations) {
0202 LOGDRESSED("Basic2DGenericPFlowClusterizer:growAndStabilizePFClusters")
0203 << "reached " << _maxIterations << " iterations, terminated position "
0204 << "fit with diff = " << diff;
0205 }
0206 if (iter >= _maxIterations || diff <= _stoppingTolerance * toleranceScaling)
0207 return;
0208
0209 std::vector<reco::PFCluster::REPPoint> clus_prev_pos;
0210 for (auto& cluster : clusters) {
0211 const reco::PFCluster::REPPoint& repp = cluster.positionREP();
0212 clus_prev_pos.emplace_back(repp.rho(), repp.eta(), repp.phi());
0213 if (_convergencePosCalc) {
0214 if (clusters.size() == 1 && _allCellsPosCalc) {
0215 _allCellsPosCalc->calculateAndSetPosition(cluster, hcalCuts);
0216 } else {
0217 _positionCalc->calculateAndSetPosition(cluster, hcalCuts);
0218 }
0219 }
0220 cluster.resetHitsAndFractions();
0221 }
0222
0223 std::vector<double> dist2, frac;
0224 double fractot = 0;
0225 for (const reco::PFRecHitFraction& rhf : topo.recHitFractions()) {
0226 const reco::PFRecHitRef& refhit = rhf.recHitRef();
0227 int cell_layer = (int)refhit->layer();
0228 if (cell_layer == PFLayer::HCAL_BARREL2 && std::abs(refhit->positionREP().eta()) > 0.34) {
0229 cell_layer *= 100;
0230 }
0231
0232 math::XYZPoint topocellpos_xyz(refhit->position());
0233 dist2.clear();
0234 frac.clear();
0235 fractot = 0;
0236
0237 double recHitEnergyNorm = 0.;
0238 auto const& recHitEnergyNormDepthPair = _recHitEnergyNorms.find(cell_layer)->second;
0239
0240 if (hcalCuts != nullptr &&
0241 (cell_layer == PFLayer::HCAL_BARREL1 || cell_layer == PFLayer::HCAL_ENDCAP)) {
0242 HcalDetId thisId = refhit->detId();
0243 const HcalPFCut* item = hcalCuts->getValues(thisId.rawId());
0244 recHitEnergyNorm = item->noiseThreshold();
0245 } else {
0246 for (unsigned int j = 0; j < recHitEnergyNormDepthPair.second.size(); ++j) {
0247 int depth = recHitEnergyNormDepthPair.first[j];
0248 if ((cell_layer == PFLayer::HCAL_BARREL1 && refhit->depth() == depth) ||
0249 (cell_layer == PFLayer::HCAL_ENDCAP && refhit->depth() == depth) ||
0250 (cell_layer != PFLayer::HCAL_ENDCAP && cell_layer != PFLayer::HCAL_BARREL1))
0251 recHitEnergyNorm = recHitEnergyNormDepthPair.second[j];
0252 }
0253 }
0254
0255
0256 for (auto& cluster : clusters) {
0257 const math::XYZPoint& clusterpos_xyz = cluster.position();
0258 const math::XYZVector deltav = clusterpos_xyz - topocellpos_xyz;
0259 const double d2 = deltav.Mag2() / _showerSigma2;
0260 dist2.emplace_back(d2);
0261 if (d2 > 100) {
0262 LOGDRESSED("Basic2DGenericPFlowClusterizer:growAndStabilizePFClusters")
0263 << "Warning! :: pfcluster-topocell distance is too large! d= " << d2;
0264 }
0265
0266
0267 double fraction;
0268 if (refhit->detId() == cluster.seed() && _excludeOtherSeeds) {
0269 fraction = 1.0;
0270 } else if (seedable[refhit.key()] && _excludeOtherSeeds) {
0271 fraction = 0.0;
0272 } else {
0273 fraction = cluster.energy() / recHitEnergyNorm * vdt::fast_expf(-0.5 * d2);
0274 }
0275 fractot += fraction;
0276 frac.emplace_back(fraction);
0277 }
0278 for (unsigned i = 0; i < clusters.size(); ++i) {
0279 if (fractot > _minFracTot || (refhit->detId() == clusters[i].seed() && fractot > 0.0)) {
0280 frac[i] /= fractot;
0281 } else {
0282 continue;
0283 }
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295 if (dist2[i] < 100.0 || frac[i] > 0.9999) {
0296 clusters[i].addRecHitFraction(reco::PFRecHitFraction(refhit, frac[i]));
0297 }
0298 }
0299 }
0300
0301 double diff2 = 0.0;
0302 for (unsigned i = 0; i < clusters.size(); ++i) {
0303 if (_convergencePosCalc) {
0304 _convergencePosCalc->calculateAndSetPosition(clusters[i], hcalCuts);
0305 } else {
0306 if (clusters.size() == 1 && _allCellsPosCalc) {
0307 _allCellsPosCalc->calculateAndSetPosition(clusters[i], hcalCuts);
0308 } else {
0309 _positionCalc->calculateAndSetPosition(clusters[i], hcalCuts);
0310 }
0311 }
0312 const double delta2 = reco::deltaR2(clusters[i].positionREP(), clus_prev_pos[i]);
0313 if (delta2 > diff2)
0314 diff2 = delta2;
0315 }
0316 diff = std::sqrt(diff2);
0317 dist2.clear();
0318 frac.clear();
0319 clus_prev_pos.clear();
0320 growPFClusters(topo, seedable, toleranceScaling, iter + 1, diff, clusters, hcalCuts);
0321 }
0322
0323 void Basic2DGenericPFlowClusterizer::prunePFClusters(reco::PFClusterCollection& clusters) const {
0324 for (auto& cluster : clusters) {
0325 cluster.pruneUsing([&](const reco::PFRecHitFraction& rhf) { return rhf.fraction() > _minFractionToKeep; });
0326 }
0327 }