PFClusterProducer/plugins/Cluster3DPCACalculator.cc

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) {  // time == -1 means no measurement
0073       // all times are offset by one nanosecond in digitizer
0074       // remove that here so all times of flight
0075       // are with respect to (0,0,0)
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       //temporarily changed exception to warning
0089       //      throw cms::Exception("PFClusterAlgo")
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);  // put rec_hit energy in units of 10 MeV
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   // calculate the position
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 }