PFClusterProducer/plugins/PFlow2DClusterizerWithTime.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/CaloRecHitResolutionProvider.h"
0007 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterBuilderBase.h"
0008 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0009 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0010
0011 #include "Math/GenVector/VectorUtil.h"
0012 #include "TMath.h"
0013 #include "vdt/vdtMath.h"
0014
0015 #include <iterator>
0016 #include <unordered_map>
0017
0018 class PFlow2DClusterizerWithTime : public PFClusterBuilderBase {
0019   typedef PFlow2DClusterizerWithTime B2DGPF;
0020
0021 public:
0022   PFlow2DClusterizerWithTime(const edm::ParameterSet& conf, edm::ConsumesCollector& cc);
0023
0024   ~PFlow2DClusterizerWithTime() override = default;
0025   PFlow2DClusterizerWithTime(const B2DGPF&) = delete;
0026   B2DGPF& operator=(const B2DGPF&) = delete;
0027
0028   void update(const edm::EventSetup& es) override {
0029     _positionCalc->update(es);
0030     if (_allCellsPosCalc)
0031       _allCellsPosCalc->update(es);
0032     if (_convergencePosCalc)
0033       _convergencePosCalc->update(es);
0034   }
0035
0036   void buildClusters(const reco::PFClusterCollection&,
0037                      const std::vector<bool>&,
0038                      reco::PFClusterCollection& outclus,
0039                      const HcalPFCuts*) override;
0040
0041 private:
0042   const unsigned _maxIterations;
0043   const double _stoppingTolerance;
0044   const double _showerSigma2;
0045   const double _timeSigma_eb;
0046   const double _timeSigma_ee;
0047   const bool _excludeOtherSeeds;
0048   const double _minFracTot;
0049   const double _maxNSigmaTime;
0050   const double _minChi2Prob;
0051   const bool _clusterTimeResFromSeed;
0052
0053   const std::unordered_map<std::string, int> _layerMap;
0054   std::unordered_map<int, double> _recHitEnergyNorms;
0055   std::unique_ptr<PFCPositionCalculatorBase> _allCellsPosCalc;
0056   std::unique_ptr<PFCPositionCalculatorBase> _convergencePosCalc;
0057
0058   std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcBarrel;
0059   std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcEndcap;
0060
0061   void seedPFClustersFromTopo(const reco::PFCluster&,
0062                               const std::vector<bool>&,
0063                               reco::PFClusterCollection&,
0064                               const HcalPFCuts*) const;
0065
0066   void growPFClusters(const reco::PFCluster&,
0067                       const std::vector<bool>&,
0068                       const unsigned toleranceScaling,
0069                       const unsigned iter,
0070                       double dist,
0071                       reco::PFClusterCollection&,
0072                       const HcalPFCuts* hcalCuts) const;
0073   void clusterTimeResolution(reco::PFCluster& cluster, double& res) const;
0074   void clusterTimeResolutionFromSeed(reco::PFCluster& cluster, double& res) const;
0075   double dist2Time(const reco::PFCluster&, const reco::PFRecHitRef&, int cell_layer, double prev_timeres2) const;
0076   bool passChi2Prob(size_t iCluster,
0077                     double dist2,
0078                     std::vector<double>& clus_chi2,
0079                     std::vector<size_t>& clus_chi2_nhits) const;
0080   void prunePFClusters(reco::PFClusterCollection&) const;
0081 };
0082
0083 DEFINE_EDM_PLUGIN(PFClusterBuilderFactory, PFlow2DClusterizerWithTime, "PFlow2DClusterizerWithTime");
0084
0085 #ifdef PFLOW_DEBUG
0086 #define LOGVERB(x) edm::LogVerbatim(x)
0087 #define LOGWARN(x) edm::LogWarning(x)
0088 #define LOGERR(x) edm::LogError(x)
0089 #define LOGDRESSED(x) edm::LogInfo(x)
0090 #else
0091 #define LOGVERB(x) LogTrace(x)
0092 #define LOGWARN(x) edm::LogWarning(x)
0093 #define LOGERR(x) edm::LogError(x)
0094 #define LOGDRESSED(x) LogDebug(x)
0095 #endif
0096
0097 PFlow2DClusterizerWithTime::PFlow2DClusterizerWithTime(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0098     : PFClusterBuilderBase(conf, cc),
0099       _maxIterations(conf.getParameter<unsigned>("maxIterations")),
0100       _stoppingTolerance(conf.getParameter<double>("stoppingTolerance")),
0101       _showerSigma2(std::pow(conf.getParameter<double>("showerSigma"), 2.0)),
0102       _timeSigma_eb(std::pow(conf.getParameter<double>("timeSigmaEB"), 2.0)),
0103       _timeSigma_ee(std::pow(conf.getParameter<double>("timeSigmaEE"), 2.0)),
0104       _excludeOtherSeeds(conf.getParameter<bool>("excludeOtherSeeds")),
0105       _minFracTot(conf.getParameter<double>("minFracTot")),
0106       _maxNSigmaTime(std::pow(conf.getParameter<double>("maxNSigmaTime"), 2.0)),
0107       _minChi2Prob(conf.getParameter<double>("minChi2Prob")),
0108       _clusterTimeResFromSeed(conf.getParameter<bool>("clusterTimeResFromSeed")),
0109
0110       _layerMap({{"PS2", (int)PFLayer::PS2},
0111                  {"PS1", (int)PFLayer::PS1},
0112                  {"ECAL_ENDCAP", (int)PFLayer::ECAL_ENDCAP},
0113                  {"ECAL_BARREL", (int)PFLayer::ECAL_BARREL},
0114                  {"NONE", (int)PFLayer::NONE},
0115                  {"HCAL_BARREL1", (int)PFLayer::HCAL_BARREL1},
0116                  {"HCAL_BARREL2_RING0", (int)PFLayer::HCAL_BARREL2},
0117                  {"HCAL_BARREL2_RING1", 100 * (int)PFLayer::HCAL_BARREL2},
0118                  {"HCAL_ENDCAP", (int)PFLayer::HCAL_ENDCAP},
0119                  {"HF_EM", (int)PFLayer::HF_EM},
0120                  {"HF_HAD", (int)PFLayer::HF_HAD}}) {
0121   const std::vector<edm::ParameterSet>& thresholds = conf.getParameterSetVector("recHitEnergyNorms");
0122   for (const auto& pset : thresholds) {
0123     const std::string& det = pset.getParameter<std::string>("detector");
0124     const double& rhE_norm = pset.getParameter<double>("recHitEnergyNorm");
0125     auto entry = _layerMap.find(det);
0126     if (entry == _layerMap.end()) {
0127       throw cms::Exception("InvalidDetectorLayer") << "Detector layer : " << det << " is not in the list of recognized"
0128                                                    << " detector layers!";
0129     }
0130     _recHitEnergyNorms.emplace(_layerMap.find(det)->second, rhE_norm);
0131   }
0132
0133   const auto& acConf = conf.getParameterSet("allCellsPositionCalc");
0134   if (!acConf.empty()) {
0135     const std::string& algoac = acConf.getParameter<std::string>("algoName");
0136     if (!algoac.empty())
0137       _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0138   }
0139   // if necessary a third pos calc for convergence testing
0140   const auto& convConf = conf.getParameterSet("positionCalcForConvergence");
0141   if (!convConf.empty()) {
0142     const std::string& algoconv = convConf.getParameter<std::string>("algoName");
0143     if (!algoconv.empty())
0144       _convergencePosCalc = PFCPositionCalculatorFactory::get()->create(algoconv, convConf, cc);
0145   }
0146   const auto& timeResConfBarrel = conf.getParameterSet("timeResolutionCalcBarrel");
0147   if (!timeResConfBarrel.empty())
0148     _timeResolutionCalcBarrel = std::make_unique<CaloRecHitResolutionProvider>(timeResConfBarrel);
0149   const auto& timeResConfEndcap = conf.getParameterSet("timeResolutionCalcEndcap");
0150   if (!timeResConfEndcap.empty())
0151     _timeResolutionCalcEndcap = std::make_unique<CaloRecHitResolutionProvider>(timeResConfEndcap);
0152 }
0153
0154 void PFlow2DClusterizerWithTime::buildClusters(const reco::PFClusterCollection& input,
0155                                                const std::vector<bool>& seedable,
0156                                                reco::PFClusterCollection& output,
0157                                                const HcalPFCuts* hcalCuts) {
0158   reco::PFClusterCollection clustersInTopo;
0159   for (const auto& topocluster : input) {
0160     clustersInTopo.clear();
0161     seedPFClustersFromTopo(topocluster, seedable, clustersInTopo, hcalCuts);
0162     const unsigned tolScal = std::pow(std::max(1.0, clustersInTopo.size() - 1.0), 2.0);
0163     growPFClusters(topocluster, seedable, tolScal, 0, tolScal, clustersInTopo, hcalCuts);
0164     // step added by Josh Bendavid, removes low-fraction clusters
0165     // did not impact position resolution with fraction cut of 1e-7
0166     // decreases the size of each pf cluster considerably
0167     prunePFClusters(clustersInTopo);
0168     // recalculate the positions of the pruned clusters
0169     if (_convergencePosCalc) {
0170       // if defined, use the special position calculation for convergence tests
0171       _convergencePosCalc->calculateAndSetPositions(clustersInTopo, hcalCuts);
0172     } else {
0173       if (clustersInTopo.size() == 1 && _allCellsPosCalc) {
0174         _allCellsPosCalc->calculateAndSetPosition(clustersInTopo.back(), hcalCuts);
0175       } else {
0176         _positionCalc->calculateAndSetPositions(clustersInTopo, hcalCuts);
0177       }
0178     }
0179     for (auto& clusterout : clustersInTopo) {
0180       output.insert(output.end(), std::move(clusterout));
0181     }
0182   }
0183 }
0184
0185 void PFlow2DClusterizerWithTime::seedPFClustersFromTopo(const reco::PFCluster& topo,
0186                                                         const std::vector<bool>& seedable,
0187                                                         reco::PFClusterCollection& initialPFClusters,
0188                                                         const HcalPFCuts* hcalCuts) const {
0189   const auto& recHitFractions = topo.recHitFractions();
0190   for (const auto& rhf : recHitFractions) {
0191     if (!seedable[rhf.recHitRef().key()])
0192       continue;
0193     initialPFClusters.push_back(reco::PFCluster());
0194     reco::PFCluster& current = initialPFClusters.back();
0195     current.addRecHitFraction(rhf);
0196     current.setSeed(rhf.recHitRef()->detId());
0197     if (_convergencePosCalc) {
0198       _convergencePosCalc->calculateAndSetPosition(current, hcalCuts);
0199     } else {
0200       _positionCalc->calculateAndSetPosition(current, hcalCuts);
0201     }
0202   }
0203 }
0204
0205 void PFlow2DClusterizerWithTime::growPFClusters(const reco::PFCluster& topo,
0206                                                 const std::vector<bool>& seedable,
0207                                                 const unsigned toleranceScaling,
0208                                                 const unsigned iter,
0209                                                 double diff,
0210                                                 reco::PFClusterCollection& clusters,
0211                                                 const HcalPFCuts* hcalCuts) const {
0212   if (iter >= _maxIterations) {
0213     LOGDRESSED("PFlow2DClusterizerWithTime:growAndStabilizePFClusters")
0214         << "reached " << _maxIterations << " iterations, terminated position "
0215         << "fit with diff = " << diff;
0216   }
0217   if (iter >= _maxIterations || diff <= _stoppingTolerance * toleranceScaling)
0218     return;
0219   // reset the rechits in this cluster, keeping the previous position
0220   std::vector<reco::PFCluster::REPPoint> clus_prev_pos;
0221   // also calculate and keep the previous time resolution
0222   std::vector<double> clus_prev_timeres2;
0223
0224   for (auto& cluster : clusters) {
0225     const reco::PFCluster::REPPoint& repp = cluster.positionREP();
0226     clus_prev_pos.emplace_back(repp.rho(), repp.eta(), repp.phi());
0227     if (_convergencePosCalc) {
0228       if (clusters.size() == 1 && _allCellsPosCalc) {
0229         _allCellsPosCalc->calculateAndSetPosition(cluster, hcalCuts);
0230       } else {
0231         _positionCalc->calculateAndSetPosition(cluster, hcalCuts);
0232       }
0233     }
0234     double resCluster2;
0235     if (_clusterTimeResFromSeed)
0236       clusterTimeResolutionFromSeed(cluster, resCluster2);
0237     else
0238       clusterTimeResolution(cluster, resCluster2);
0239     clus_prev_timeres2.push_back(resCluster2);
0240     cluster.resetHitsAndFractions();
0241   }
0242   // loop over topo cluster and grow current PFCluster hypothesis
0243   std::vector<double> dist2, frac;
0244
0245   // Store chi2 values and nhits to calculate chi2 prob. inline
0246   std::vector<double> clus_chi2;
0247   std::vector<size_t> clus_chi2_nhits;
0248
0249   double fractot = 0;
0250   for (const reco::PFRecHitFraction& rhf : topo.recHitFractions()) {
0251     const reco::PFRecHitRef& refhit = rhf.recHitRef();
0252     int cell_layer = (int)refhit->layer();
0253     if (cell_layer == PFLayer::HCAL_BARREL2 && std::abs(refhit->positionREP().eta()) > 0.34) {
0254       cell_layer *= 100;
0255     }
0256     double recHitEnergyNorm = _recHitEnergyNorms.find(cell_layer)->second;
0257     if (hcalCuts != nullptr) {  // this means, cutsFromDB is set to True in the producer code
0258       if ((cell_layer == PFLayer::HCAL_BARREL1) || (cell_layer == PFLayer::HCAL_ENDCAP)) {
0259         HcalDetId thisId = refhit->detId();
0260         const HcalPFCut* item = hcalCuts->getValues(thisId.rawId());
0261         recHitEnergyNorm = item->noiseThreshold();
0262       }
0263     }
0264
0265     math::XYZPoint topocellpos_xyz(refhit->position());
0266     dist2.clear();
0267     frac.clear();
0268     fractot = 0;
0269
0270     // add rechits to clusters, calculating fraction based on distance
0271     // need position in vector to get cluster time resolution
0272     for (size_t iCluster = 0; iCluster < clusters.size(); ++iCluster) {
0273       reco::PFCluster& cluster = clusters[iCluster];
0274       const math::XYZPoint& clusterpos_xyz = cluster.position();
0275       const math::XYZVector deltav = clusterpos_xyz - topocellpos_xyz;
0276       double d2 = deltav.Mag2() / _showerSigma2;
0277
0278       double d2time = dist2Time(cluster, refhit, cell_layer, clus_prev_timeres2[iCluster]);
0279       d2 += d2time;
0280
0281       if (_minChi2Prob > 0. && !passChi2Prob(iCluster, d2time, clus_chi2, clus_chi2_nhits))
0282         d2 = 999.;
0283
0284       dist2.emplace_back(d2);
0285
0286       if (d2 > 100) {
0287         LOGDRESSED("PFlow2DClusterizerWithTime:growAndStabilizePFClusters")
0288             << "Warning! :: pfcluster-topocell distance is too large! d= " << d2;
0289       }
0290       // fraction assignment logic
0291       double fraction;
0292       if (refhit->detId() == cluster.seed() && _excludeOtherSeeds) {
0293         fraction = 1.0;
0294       } else if (seedable[refhit.key()] && _excludeOtherSeeds) {
0295         fraction = 0.0;
0296       } else {
0297         fraction = cluster.energy() / recHitEnergyNorm * vdt::fast_expf(-0.5 * d2);
0298       }
0299       fractot += fraction;
0300       frac.emplace_back(fraction);
0301     }
0302     for (unsigned i = 0; i < clusters.size(); ++i) {
0303       if (fractot > _minFracTot || (refhit->detId() == clusters[i].seed() && fractot > 0.0)) {
0304         frac[i] /= fractot;
0305       } else {
0306         continue;
0307       }
0308       // if the fraction has been set to 0, the cell
0309       // is now added to the cluster - careful ! (PJ, 19/07/08)
0310       // BUT KEEP ONLY CLOSE CELLS OTHERWISE MEMORY JUST EXPLOSES
0311       // (PJ, 15/09/08 <- similar to what existed before the
0312       // previous bug fix, but keeps the close seeds inside,
0313       // even if their fraction was set to zero.)
0314       // Also add a protection to keep the seed in the cluster
0315       // when the latter gets far from the former. These cases
0316       // (about 1% of the clusters) need to be studied, as
0317       // they create fake photons, in general.
0318       // (PJ, 16/09/08)
0319       if (dist2[i] < 100.0 || frac[i] > 0.9999) {
0320         clusters[i].addRecHitFraction(reco::PFRecHitFraction(refhit, frac[i]));
0321       }
0322     }
0323   }
0324   // recalculate positions and calculate convergence parameter
0325   double diff2 = 0.0;
0326   for (unsigned i = 0; i < clusters.size(); ++i) {
0327     if (_convergencePosCalc) {
0328       _convergencePosCalc->calculateAndSetPosition(clusters[i], hcalCuts);
0329     } else {
0330       if (clusters.size() == 1 && _allCellsPosCalc) {
0331         _allCellsPosCalc->calculateAndSetPosition(clusters[i], hcalCuts);
0332       } else {
0333         _positionCalc->calculateAndSetPosition(clusters[i], hcalCuts);
0334       }
0335     }
0336     const double delta2 = reco::deltaR2(clusters[i].positionREP(), clus_prev_pos[i]);
0337     if (delta2 > diff2)
0338       diff2 = delta2;
0339   }
0340   diff = std::sqrt(diff2);
0341   dist2.clear();
0342   frac.clear();
0343   clus_prev_pos.clear();  // avoid badness
0344   clus_chi2.clear();
0345   clus_chi2_nhits.clear();
0346   clus_prev_timeres2.clear();
0347   growPFClusters(topo, seedable, toleranceScaling, iter + 1, diff, clusters, hcalCuts);
0348 }
0349
0350 void PFlow2DClusterizerWithTime::prunePFClusters(reco::PFClusterCollection& clusters) const {
0351   for (auto& cluster : clusters) {
0352     cluster.pruneUsing([&](const reco::PFRecHitFraction& rhf) { return rhf.fraction() > _minFractionToKeep; });
0353   }
0354 }
0355
0356 void PFlow2DClusterizerWithTime::clusterTimeResolutionFromSeed(reco::PFCluster& cluster, double& clusterRes2) const {
0357   clusterRes2 = 10000.;
0358   for (auto& rhf : cluster.recHitFractions()) {
0359     const reco::PFRecHit& rh = *(rhf.recHitRef());
0360     if (rh.detId() == cluster.seed()) {
0361       cluster.setTime(rh.time());
0362       bool isBarrel = (rh.layer() == PFLayer::HCAL_BARREL1 || rh.layer() == PFLayer::HCAL_BARREL2 ||
0363                        rh.layer() == PFLayer::ECAL_BARREL);
0364       if (isBarrel)
0365         clusterRes2 = _timeResolutionCalcBarrel->timeResolution2(rh.energy());
0366       else
0367         clusterRes2 = _timeResolutionCalcEndcap->timeResolution2(rh.energy());
0368       cluster.setTimeError(std::sqrt(float(clusterRes2)));
0369     }
0370   }
0371 }
0372
0373 void PFlow2DClusterizerWithTime::clusterTimeResolution(reco::PFCluster& cluster, double& clusterRes2) const {
0374   double sumTimeSigma2 = 0.;
0375   double sumSigma2 = 0.;
0376
0377   for (auto& rhf : cluster.recHitFractions()) {
0378     const reco::PFRecHit& rh = *(rhf.recHitRef());
0379     const double rhf_f = rhf.fraction();
0380
0381     if (rhf_f == 0.)
0382       continue;
0383
0384     bool isBarrel = (rh.layer() == PFLayer::HCAL_BARREL1 || rh.layer() == PFLayer::HCAL_BARREL2 ||
0385                      rh.layer() == PFLayer::ECAL_BARREL);
0386     double res2 = 10000.;
0387     if (isBarrel)
0388       res2 = _timeResolutionCalcBarrel->timeResolution2(rh.energy());
0389     else
0390       res2 = _timeResolutionCalcEndcap->timeResolution2(rh.energy());
0391
0392     sumTimeSigma2 += rhf_f * rh.time() / res2;
0393     sumSigma2 += rhf_f / res2;
0394   }
0395   if (sumSigma2 > 0.) {
0396     clusterRes2 = 1. / sumSigma2;
0397     cluster.setTime(sumTimeSigma2 / sumSigma2);
0398     cluster.setTimeError(std::sqrt(float(clusterRes2)));
0399   } else {
0400     clusterRes2 = 999999.;
0401   }
0402 }
0403
0404 double PFlow2DClusterizerWithTime::dist2Time(const reco::PFCluster& cluster,
0405                                              const reco::PFRecHitRef& refhit,
0406                                              int cell_layer,
0407                                              double prev_timeres2) const {
0408   const double deltaT = cluster.time() - refhit->time();
0409   const double t2 = deltaT * deltaT;
0410   double res2 = 100.;
0411
0412   if (cell_layer == PFLayer::HCAL_BARREL1 || cell_layer == PFLayer::HCAL_BARREL2 ||
0413       cell_layer == PFLayer::ECAL_BARREL) {
0414     if (_timeResolutionCalcBarrel) {
0415       const double resCluster2 = prev_timeres2;
0416       res2 = resCluster2 + _timeResolutionCalcBarrel->timeResolution2(refhit->energy());
0417     } else {
0418       return t2 / _timeSigma_eb;
0419     }
0420   } else if (cell_layer == PFLayer::HCAL_ENDCAP || cell_layer == PFLayer::HF_EM || cell_layer == PFLayer::HF_HAD ||
0421              cell_layer == PFLayer::ECAL_ENDCAP) {
0422     if (_timeResolutionCalcEndcap) {
0423       const double resCluster2 = prev_timeres2;
0424       res2 = resCluster2 + _timeResolutionCalcEndcap->timeResolution2(refhit->energy());
0425     } else {
0426       return t2 / _timeSigma_ee;
0427     }
0428   }
0429
0430   double distTime2 = t2 / res2;
0431   if (distTime2 > _maxNSigmaTime)
0432     return 999.;  // reject hit
0433
0434   return distTime2;
0435 }
0436
0437 bool PFlow2DClusterizerWithTime::passChi2Prob(size_t iCluster,
0438                                               double dist2,
0439                                               std::vector<double>& clus_chi2,
0440                                               std::vector<size_t>& clus_chi2_nhits) const {
0441   if (iCluster >= clus_chi2.size()) {  // first hit
0442     clus_chi2.push_back(dist2);
0443     clus_chi2_nhits.push_back(1);
0444     return true;
0445   } else {
0446     double chi2 = clus_chi2[iCluster];
0447     size_t nhitsCluster = clus_chi2_nhits[iCluster];
0448     chi2 += dist2;
0449     if (TMath::Prob(chi2, nhitsCluster) >= _minChi2Prob) {
0450       clus_chi2[iCluster] = chi2;
0451       clus_chi2_nhits[iCluster] = nhitsCluster + 1;
0452       return true;
0453     }
0454   }
0455   return false;
0456 }