Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // if necessary a third pos calc for convergence testing
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     // step added by Josh Bendavid, removes low-fraction clusters
0157     // did not impact position resolution with fraction cut of 1e-7
0158     // decreases the size of each pf cluster considerably
0159     prunePFClusters(clustersInTopo);
0160     // recalculate the positions of the pruned clusters
0161     if (_convergencePosCalc) {
0162       // if defined, use the special position calculation for convergence tests
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   // reset the rechits in this cluster, keeping the previous position
0210   std::vector<reco::PFCluster::REPPoint> clus_prev_pos;
0211   // also calculate and keep the previous time resolution
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   // loop over topo cluster and grow current PFCluster hypothesis
0233   std::vector<double> dist2, frac;
0234 
0235   // Store chi2 values and nhits to calculate chi2 prob. inline
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     // add rechits to clusters, calculating fraction based on distance
0253     // need position in vector to get cluster time resolution
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       // fraction assignment logic
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       // if the fraction has been set to 0, the cell
0291       // is now added to the cluster - careful ! (PJ, 19/07/08)
0292       // BUT KEEP ONLY CLOSE CELLS OTHERWISE MEMORY JUST EXPLOSES
0293       // (PJ, 15/09/08 <- similar to what existed before the
0294       // previous bug fix, but keeps the close seeds inside,
0295       // even if their fraction was set to zero.)
0296       // Also add a protection to keep the seed in the cluster
0297       // when the latter gets far from the former. These cases
0298       // (about 1% of the clusters) need to be studied, as
0299       // they create fake photons, in general.
0300       // (PJ, 16/09/08)
0301       if (dist2[i] < 100.0 || frac[i] > 0.9999) {
0302         clusters[i].addRecHitFraction(reco::PFRecHitFraction(refhit, frac[i]));
0303       }
0304     }
0305   }
0306   // recalculate positions and calculate convergence parameter
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();  // avoid badness
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.;  // reject hit
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()) {  // first hit
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!