File indexing completed on 2023-03-17 11:21:03
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003
0004 #include <cmath>
0005 #include <memory>
0006
0007 #include <unordered_map>
0008
0009 #include "vdt/vdtMath.h"
0010
0011 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0012 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0013 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015
0016 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0017
0018 #include "TPrincipal.h"
0019
0020 class Cluster3DPCACalculator : public PFCPositionCalculatorBase {
0021 public:
0022 Cluster3DPCACalculator(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0023 : PFCPositionCalculatorBase(conf, cc),
0024 updateTiming_(conf.getParameter<bool>("updateTiming")),
0025 pca_(new TPrincipal(3, "D")) {}
0026 Cluster3DPCACalculator(const Cluster3DPCACalculator&) = delete;
0027 Cluster3DPCACalculator& operator=(const Cluster3DPCACalculator&) = delete;
0028
0029 void calculateAndSetPosition(reco::PFCluster&) override;
0030 void calculateAndSetPositions(reco::PFClusterCollection&) override;
0031
0032 private:
0033 const bool updateTiming_;
0034 std::unique_ptr<TPrincipal> pca_;
0035
0036 void showerParameters(const reco::PFCluster&, math::XYZPoint&, math::XYZVector&);
0037
0038 void calculateAndSetPositionActual(reco::PFCluster&);
0039 };
0040
0041 DEFINE_EDM_PLUGIN(PFCPositionCalculatorFactory, Cluster3DPCACalculator, "Cluster3DPCACalculator");
0042
0043 void Cluster3DPCACalculator::calculateAndSetPosition(reco::PFCluster& cluster) {
0044 pca_ = std::make_unique<TPrincipal>(3, "D");
0045 calculateAndSetPositionActual(cluster);
0046 }
0047
0048 void Cluster3DPCACalculator::calculateAndSetPositions(reco::PFClusterCollection& clusters) {
0049 for (reco::PFCluster& cluster : clusters) {
0050 pca_ = std::make_unique<TPrincipal>(3, "D");
0051 calculateAndSetPositionActual(cluster);
0052 }
0053 }
0054
0055 void Cluster3DPCACalculator::calculateAndSetPositionActual(reco::PFCluster& cluster) {
0056 if (!cluster.seed()) {
0057 throw cms::Exception("ClusterWithNoSeed") << " Found a cluster with no seed: " << cluster;
0058 }
0059 double cl_energy = 0;
0060 double max_e = 0.0;
0061 double avg_time = 0.0;
0062 double time_norm = 0.0;
0063 PFLayer::Layer max_e_layer = PFLayer::NONE;
0064 reco::PFRecHitRef refseed;
0065 double pcavars[3];
0066
0067 for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
0068 const reco::PFRecHitRef& refhit = rhf.recHitRef();
0069 double rh_energy = refhit->energy();
0070 double rh_time = refhit->time();
0071 cl_energy += rh_energy * rhf.fraction();
0072 if (rh_time > 0.0) {
0073
0074
0075
0076 avg_time += (rh_time - 1.0);
0077 time_norm += 1.0;
0078 }
0079 if (rh_energy > max_e) {
0080 max_e = rh_energy;
0081 max_e_layer = rhf.recHitRef()->layer();
0082 }
0083 if (refhit->detId() == cluster.seed())
0084 refseed = refhit;
0085 const double rh_fraction = rhf.fraction();
0086 rh_energy = refhit->energy() * rh_fraction;
0087 if (edm::isNotFinite(rh_energy)) {
0088
0089
0090 edm::LogWarning("PFClusterAlgo") << "rechit " << refhit->detId() << " has a NaN energy... "
0091 << "The input of the particle flow clustering seems to be corrupted.";
0092 continue;
0093 }
0094 pcavars[0] = refhit->position().x();
0095 pcavars[1] = refhit->position().y();
0096 pcavars[2] = refhit->position().z();
0097 int nhit = int(rh_energy * 100);
0098
0099 for (int i = 0; i < nhit; ++i) {
0100 pca_->AddRow(pcavars);
0101 }
0102 }
0103 cluster.setEnergy(cl_energy);
0104 cluster.setLayer(max_e_layer);
0105
0106
0107 pca_->MakePrincipals();
0108 const TVectorD& means = *(pca_->GetMeanValues());
0109 const TMatrixD& eigens = *(pca_->GetEigenVectors());
0110
0111 math::XYZPoint barycenter(means[0], means[1], means[2]);
0112 math::XYZVector axis(eigens(0, 0), eigens(1, 0), eigens(2, 0));
0113
0114 if (time_norm > 0.0) {
0115 avg_time = avg_time / time_norm;
0116 } else {
0117 avg_time = std::numeric_limits<double>::min();
0118 }
0119
0120 if (axis.z() * barycenter.z() < 0.0) {
0121 axis = math::XYZVector(-eigens(0, 0), -eigens(1, 0), -eigens(2, 0));
0122 }
0123
0124 if (updateTiming_) {
0125 cluster.setTime(avg_time);
0126 }
0127 cluster.setPosition(barycenter);
0128 cluster.calculatePositionREP();
0129 }