File indexing completed on 2021-07-30 02:33:27
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/CaloRecHitResolutionProvider.h"
0007 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterBuilderBase.h"
0008
0009 #include "Math/GenVector/VectorUtil.h"
0010 #include "TMath.h"
0011 #include "vdt/vdtMath.h"
0012
0013 #include <iterator>
0014 #include <unordered_map>
0015
0016 class PFlow2DClusterizerWithTime : public PFClusterBuilderBase {
0017 typedef PFlow2DClusterizerWithTime B2DGPF;
0018
0019 public:
0020 PFlow2DClusterizerWithTime(const edm::ParameterSet& conf, edm::ConsumesCollector& cc);
0021
0022 ~PFlow2DClusterizerWithTime() override = default;
0023 PFlow2DClusterizerWithTime(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) override;
0037
0038 private:
0039 const unsigned _maxIterations;
0040 const double _stoppingTolerance;
0041 const double _showerSigma2;
0042 const double _timeSigma_eb;
0043 const double _timeSigma_ee;
0044 const bool _excludeOtherSeeds;
0045 const double _minFracTot;
0046 const double _maxNSigmaTime;
0047 const double _minChi2Prob;
0048 const bool _clusterTimeResFromSeed;
0049
0050 const std::unordered_map<std::string, int> _layerMap;
0051 std::unordered_map<int, double> _recHitEnergyNorms;
0052 std::unique_ptr<PFCPositionCalculatorBase> _allCellsPosCalc;
0053 std::unique_ptr<PFCPositionCalculatorBase> _convergencePosCalc;
0054
0055 std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcBarrel;
0056 std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcEndcap;
0057
0058 void seedPFClustersFromTopo(const reco::PFCluster&, const std::vector<bool>&, reco::PFClusterCollection&) const;
0059
0060 void growPFClusters(const reco::PFCluster&,
0061 const std::vector<bool>&,
0062 const unsigned toleranceScaling,
0063 const unsigned iter,
0064 double dist,
0065 reco::PFClusterCollection&) const;
0066 void clusterTimeResolution(reco::PFCluster& cluster, double& res) const;
0067 void clusterTimeResolutionFromSeed(reco::PFCluster& cluster, double& res) const;
0068 double dist2Time(const reco::PFCluster&, const reco::PFRecHitRef&, int cell_layer, double prev_timeres2) const;
0069 bool passChi2Prob(size_t iCluster,
0070 double dist2,
0071 std::vector<double>& clus_chi2,
0072 std::vector<size_t>& clus_chi2_nhits) const;
0073 void prunePFClusters(reco::PFClusterCollection&) const;
0074 };
0075
0076 DEFINE_EDM_PLUGIN(PFClusterBuilderFactory, PFlow2DClusterizerWithTime, "PFlow2DClusterizerWithTime");
0077
0078 #ifdef PFLOW_DEBUG
0079 #define LOGVERB(x) edm::LogVerbatim(x)
0080 #define LOGWARN(x) edm::LogWarning(x)
0081 #define LOGERR(x) edm::LogError(x)
0082 #define LOGDRESSED(x) edm::LogInfo(x)
0083 #else
0084 #define LOGVERB(x) LogTrace(x)
0085 #define LOGWARN(x) edm::LogWarning(x)
0086 #define LOGERR(x) edm::LogError(x)
0087 #define LOGDRESSED(x) LogDebug(x)
0088 #endif
0089
0090 PFlow2DClusterizerWithTime::PFlow2DClusterizerWithTime(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0091 : PFClusterBuilderBase(conf, cc),
0092 _maxIterations(conf.getParameter<unsigned>("maxIterations")),
0093 _stoppingTolerance(conf.getParameter<double>("stoppingTolerance")),
0094 _showerSigma2(std::pow(conf.getParameter<double>("showerSigma"), 2.0)),
0095 _timeSigma_eb(std::pow(conf.getParameter<double>("timeSigmaEB"), 2.0)),
0096 _timeSigma_ee(std::pow(conf.getParameter<double>("timeSigmaEE"), 2.0)),
0097 _excludeOtherSeeds(conf.getParameter<bool>("excludeOtherSeeds")),
0098 _minFracTot(conf.getParameter<double>("minFracTot")),
0099 _maxNSigmaTime(std::pow(conf.getParameter<double>("maxNSigmaTime"), 2.0)),
0100 _minChi2Prob(conf.getParameter<double>("minChi2Prob")),
0101 _clusterTimeResFromSeed(conf.getParameter<bool>("clusterTimeResFromSeed")),
0102
0103 _layerMap({{"PS2", (int)PFLayer::PS2},
0104 {"PS1", (int)PFLayer::PS1},
0105 {"ECAL_ENDCAP", (int)PFLayer::ECAL_ENDCAP},
0106 {"ECAL_BARREL", (int)PFLayer::ECAL_BARREL},
0107 {"NONE", (int)PFLayer::NONE},
0108 {"HCAL_BARREL1", (int)PFLayer::HCAL_BARREL1},
0109 {"HCAL_BARREL2_RING0", (int)PFLayer::HCAL_BARREL2},
0110 {"HCAL_BARREL2_RING1", 100 * (int)PFLayer::HCAL_BARREL2},
0111 {"HCAL_ENDCAP", (int)PFLayer::HCAL_ENDCAP},
0112 {"HF_EM", (int)PFLayer::HF_EM},
0113 {"HF_HAD", (int)PFLayer::HF_HAD}}) {
0114 const std::vector<edm::ParameterSet>& thresholds = conf.getParameterSetVector("recHitEnergyNorms");
0115 for (const auto& pset : thresholds) {
0116 const std::string& det = pset.getParameter<std::string>("detector");
0117 const double& rhE_norm = pset.getParameter<double>("recHitEnergyNorm");
0118 auto entry = _layerMap.find(det);
0119 if (entry == _layerMap.end()) {
0120 throw cms::Exception("InvalidDetectorLayer") << "Detector layer : " << det << " is not in the list of recognized"
0121 << " detector layers!";
0122 }
0123 _recHitEnergyNorms.emplace(_layerMap.find(det)->second, rhE_norm);
0124 }
0125
0126 if (conf.exists("allCellsPositionCalc")) {
0127 const edm::ParameterSet& acConf = conf.getParameterSet("allCellsPositionCalc");
0128 const std::string& algoac = acConf.getParameter<std::string>("algoName");
0129 _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0130 }
0131
0132 if (conf.exists("positionCalcForConvergence")) {
0133 const edm::ParameterSet& convConf = conf.getParameterSet("positionCalcForConvergence");
0134 const std::string& algoconv = convConf.getParameter<std::string>("algoName");
0135 _convergencePosCalc = PFCPositionCalculatorFactory::get()->create(algoconv, convConf, cc);
0136 }
0137 if (conf.exists("timeResolutionCalcBarrel")) {
0138 const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcBarrel");
0139 _timeResolutionCalcBarrel = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
0140 }
0141 if (conf.exists("timeResolutionCalcEndcap")) {
0142 const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcEndcap");
0143 _timeResolutionCalcEndcap = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
0144 }
0145 }
0146
0147 void PFlow2DClusterizerWithTime::buildClusters(const reco::PFClusterCollection& input,
0148 const std::vector<bool>& seedable,
0149 reco::PFClusterCollection& output) {
0150 reco::PFClusterCollection clustersInTopo;
0151 for (const auto& topocluster : input) {
0152 clustersInTopo.clear();
0153 seedPFClustersFromTopo(topocluster, seedable, clustersInTopo);
0154 const unsigned tolScal = std::pow(std::max(1.0, clustersInTopo.size() - 1.0), 2.0);
0155 growPFClusters(topocluster, seedable, tolScal, 0, tolScal, clustersInTopo);
0156
0157
0158
0159 prunePFClusters(clustersInTopo);
0160
0161 if (_convergencePosCalc) {
0162
0163 _convergencePosCalc->calculateAndSetPositions(clustersInTopo);
0164 } else {
0165 if (clustersInTopo.size() == 1 && _allCellsPosCalc) {
0166 _allCellsPosCalc->calculateAndSetPosition(clustersInTopo.back());
0167 } else {
0168 _positionCalc->calculateAndSetPositions(clustersInTopo);
0169 }
0170 }
0171 for (auto& clusterout : clustersInTopo) {
0172 output.insert(output.end(), std::move(clusterout));
0173 }
0174 }
0175 }
0176
0177 void PFlow2DClusterizerWithTime::seedPFClustersFromTopo(const reco::PFCluster& topo,
0178 const std::vector<bool>& seedable,
0179 reco::PFClusterCollection& initialPFClusters) const {
0180 const auto& recHitFractions = topo.recHitFractions();
0181 for (const auto& rhf : recHitFractions) {
0182 if (!seedable[rhf.recHitRef().key()])
0183 continue;
0184 initialPFClusters.push_back(reco::PFCluster());
0185 reco::PFCluster& current = initialPFClusters.back();
0186 current.addRecHitFraction(rhf);
0187 current.setSeed(rhf.recHitRef()->detId());
0188 if (_convergencePosCalc) {
0189 _convergencePosCalc->calculateAndSetPosition(current);
0190 } else {
0191 _positionCalc->calculateAndSetPosition(current);
0192 }
0193 }
0194 }
0195
0196 void PFlow2DClusterizerWithTime::growPFClusters(const reco::PFCluster& topo,
0197 const std::vector<bool>& seedable,
0198 const unsigned toleranceScaling,
0199 const unsigned iter,
0200 double diff,
0201 reco::PFClusterCollection& clusters) const {
0202 if (iter >= _maxIterations) {
0203 LOGDRESSED("PFlow2DClusterizerWithTime:growAndStabilizePFClusters")
0204 << "reached " << _maxIterations << " iterations, terminated position "
0205 << "fit with diff = " << diff;
0206 }
0207 if (iter >= _maxIterations || diff <= _stoppingTolerance * toleranceScaling)
0208 return;
0209
0210 std::vector<reco::PFCluster::REPPoint> clus_prev_pos;
0211
0212 std::vector<double> clus_prev_timeres2;
0213
0214 for (auto& cluster : clusters) {
0215 const reco::PFCluster::REPPoint& repp = cluster.positionREP();
0216 clus_prev_pos.emplace_back(repp.rho(), repp.eta(), repp.phi());
0217 if (_convergencePosCalc) {
0218 if (clusters.size() == 1 && _allCellsPosCalc) {
0219 _allCellsPosCalc->calculateAndSetPosition(cluster);
0220 } else {
0221 _positionCalc->calculateAndSetPosition(cluster);
0222 }
0223 }
0224 double resCluster2;
0225 if (_clusterTimeResFromSeed)
0226 clusterTimeResolutionFromSeed(cluster, resCluster2);
0227 else
0228 clusterTimeResolution(cluster, resCluster2);
0229 clus_prev_timeres2.push_back(resCluster2);
0230 cluster.resetHitsAndFractions();
0231 }
0232
0233 std::vector<double> dist2, frac;
0234
0235
0236 std::vector<double> clus_chi2;
0237 std::vector<size_t> clus_chi2_nhits;
0238
0239 double fractot = 0;
0240 for (const reco::PFRecHitFraction& rhf : topo.recHitFractions()) {
0241 const reco::PFRecHitRef& refhit = rhf.recHitRef();
0242 int cell_layer = (int)refhit->layer();
0243 if (cell_layer == PFLayer::HCAL_BARREL2 && std::abs(refhit->positionREP().eta()) > 0.34) {
0244 cell_layer *= 100;
0245 }
0246 const double recHitEnergyNorm = _recHitEnergyNorms.find(cell_layer)->second;
0247 math::XYZPoint topocellpos_xyz(refhit->position());
0248 dist2.clear();
0249 frac.clear();
0250 fractot = 0;
0251
0252
0253
0254 for (size_t iCluster = 0; iCluster < clusters.size(); ++iCluster) {
0255 reco::PFCluster& cluster = clusters[iCluster];
0256 const math::XYZPoint& clusterpos_xyz = cluster.position();
0257 const math::XYZVector deltav = clusterpos_xyz - topocellpos_xyz;
0258 double d2 = deltav.Mag2() / _showerSigma2;
0259
0260 double d2time = dist2Time(cluster, refhit, cell_layer, clus_prev_timeres2[iCluster]);
0261 d2 += d2time;
0262
0263 if (_minChi2Prob > 0. && !passChi2Prob(iCluster, d2time, clus_chi2, clus_chi2_nhits))
0264 d2 = 999.;
0265
0266 dist2.emplace_back(d2);
0267
0268 if (d2 > 100) {
0269 LOGDRESSED("PFlow2DClusterizerWithTime:growAndStabilizePFClusters")
0270 << "Warning! :: pfcluster-topocell distance is too large! d= " << d2;
0271 }
0272
0273 double fraction;
0274 if (refhit->detId() == cluster.seed() && _excludeOtherSeeds) {
0275 fraction = 1.0;
0276 } else if (seedable[refhit.key()] && _excludeOtherSeeds) {
0277 fraction = 0.0;
0278 } else {
0279 fraction = cluster.energy() / recHitEnergyNorm * vdt::fast_expf(-0.5 * d2);
0280 }
0281 fractot += fraction;
0282 frac.emplace_back(fraction);
0283 }
0284 for (unsigned i = 0; i < clusters.size(); ++i) {
0285 if (fractot > _minFracTot || (refhit->detId() == clusters[i].seed() && fractot > 0.0)) {
0286 frac[i] /= fractot;
0287 } else {
0288 continue;
0289 }
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301 if (dist2[i] < 100.0 || frac[i] > 0.9999) {
0302 clusters[i].addRecHitFraction(reco::PFRecHitFraction(refhit, frac[i]));
0303 }
0304 }
0305 }
0306
0307 double diff2 = 0.0;
0308 for (unsigned i = 0; i < clusters.size(); ++i) {
0309 if (_convergencePosCalc) {
0310 _convergencePosCalc->calculateAndSetPosition(clusters[i]);
0311 } else {
0312 if (clusters.size() == 1 && _allCellsPosCalc) {
0313 _allCellsPosCalc->calculateAndSetPosition(clusters[i]);
0314 } else {
0315 _positionCalc->calculateAndSetPosition(clusters[i]);
0316 }
0317 }
0318 const double delta2 = reco::deltaR2(clusters[i].positionREP(), clus_prev_pos[i]);
0319 if (delta2 > diff2)
0320 diff2 = delta2;
0321 }
0322 diff = std::sqrt(diff2);
0323 dist2.clear();
0324 frac.clear();
0325 clus_prev_pos.clear();
0326 clus_chi2.clear();
0327 clus_chi2_nhits.clear();
0328 clus_prev_timeres2.clear();
0329 growPFClusters(topo, seedable, toleranceScaling, iter + 1, diff, clusters);
0330 }
0331
0332 void PFlow2DClusterizerWithTime::prunePFClusters(reco::PFClusterCollection& clusters) const {
0333 for (auto& cluster : clusters) {
0334 cluster.pruneUsing([&](const reco::PFRecHitFraction& rhf) { return rhf.fraction() > _minFractionToKeep; });
0335 }
0336 }
0337
0338 void PFlow2DClusterizerWithTime::clusterTimeResolutionFromSeed(reco::PFCluster& cluster, double& clusterRes2) const {
0339 clusterRes2 = 10000.;
0340 for (auto& rhf : cluster.recHitFractions()) {
0341 const reco::PFRecHit& rh = *(rhf.recHitRef());
0342 if (rh.detId() == cluster.seed()) {
0343 cluster.setTime(rh.time());
0344 bool isBarrel = (rh.layer() == PFLayer::HCAL_BARREL1 || rh.layer() == PFLayer::HCAL_BARREL2 ||
0345 rh.layer() == PFLayer::ECAL_BARREL);
0346 if (isBarrel)
0347 clusterRes2 = _timeResolutionCalcBarrel->timeResolution2(rh.energy());
0348 else
0349 clusterRes2 = _timeResolutionCalcEndcap->timeResolution2(rh.energy());
0350 cluster.setTimeError(std::sqrt(float(clusterRes2)));
0351 }
0352 }
0353 }
0354
0355 void PFlow2DClusterizerWithTime::clusterTimeResolution(reco::PFCluster& cluster, double& clusterRes2) const {
0356 double sumTimeSigma2 = 0.;
0357 double sumSigma2 = 0.;
0358
0359 for (auto& rhf : cluster.recHitFractions()) {
0360 const reco::PFRecHit& rh = *(rhf.recHitRef());
0361 const double rhf_f = rhf.fraction();
0362
0363 if (rhf_f == 0.)
0364 continue;
0365
0366 bool isBarrel = (rh.layer() == PFLayer::HCAL_BARREL1 || rh.layer() == PFLayer::HCAL_BARREL2 ||
0367 rh.layer() == PFLayer::ECAL_BARREL);
0368 double res2 = 10000.;
0369 if (isBarrel)
0370 res2 = _timeResolutionCalcBarrel->timeResolution2(rh.energy());
0371 else
0372 res2 = _timeResolutionCalcEndcap->timeResolution2(rh.energy());
0373
0374 sumTimeSigma2 += rhf_f * rh.time() / res2;
0375 sumSigma2 += rhf_f / res2;
0376 }
0377 if (sumSigma2 > 0.) {
0378 clusterRes2 = 1. / sumSigma2;
0379 cluster.setTime(sumTimeSigma2 / sumSigma2);
0380 cluster.setTimeError(std::sqrt(float(clusterRes2)));
0381 } else {
0382 clusterRes2 = 999999.;
0383 }
0384 }
0385
0386 double PFlow2DClusterizerWithTime::dist2Time(const reco::PFCluster& cluster,
0387 const reco::PFRecHitRef& refhit,
0388 int cell_layer,
0389 double prev_timeres2) const {
0390 const double deltaT = cluster.time() - refhit->time();
0391 const double t2 = deltaT * deltaT;
0392 double res2 = 100.;
0393
0394 if (cell_layer == PFLayer::HCAL_BARREL1 || cell_layer == PFLayer::HCAL_BARREL2 ||
0395 cell_layer == PFLayer::ECAL_BARREL) {
0396 if (_timeResolutionCalcBarrel) {
0397 const double resCluster2 = prev_timeres2;
0398 res2 = resCluster2 + _timeResolutionCalcBarrel->timeResolution2(refhit->energy());
0399 } else {
0400 return t2 / _timeSigma_eb;
0401 }
0402 } else if (cell_layer == PFLayer::HCAL_ENDCAP || cell_layer == PFLayer::HF_EM || cell_layer == PFLayer::HF_HAD ||
0403 cell_layer == PFLayer::ECAL_ENDCAP) {
0404 if (_timeResolutionCalcEndcap) {
0405 const double resCluster2 = prev_timeres2;
0406 res2 = resCluster2 + _timeResolutionCalcEndcap->timeResolution2(refhit->energy());
0407 } else {
0408 return t2 / _timeSigma_ee;
0409 }
0410 }
0411
0412 double distTime2 = t2 / res2;
0413 if (distTime2 > _maxNSigmaTime)
0414 return 999.;
0415
0416 return distTime2;
0417 }
0418
0419 bool PFlow2DClusterizerWithTime::passChi2Prob(size_t iCluster,
0420 double dist2,
0421 std::vector<double>& clus_chi2,
0422 std::vector<size_t>& clus_chi2_nhits) const {
0423 if (iCluster >= clus_chi2.size()) {
0424 clus_chi2.push_back(dist2);
0425 clus_chi2_nhits.push_back(1);
0426 return true;
0427 } else {
0428 double chi2 = clus_chi2[iCluster];
0429 size_t nhitsCluster = clus_chi2_nhits[iCluster];
0430 chi2 += dist2;
0431 if (TMath::Prob(chi2, nhitsCluster) >= _minChi2Prob) {
0432 clus_chi2[iCluster] = chi2;
0433 clus_chi2_nhits[iCluster] = nhitsCluster + 1;
0434 return true;
0435 }
0436 }
0437 return false;
0438 }