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 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0016 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0017 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0018
0019 #include "TPrincipal.h"
0020
0021 class Cluster3DPCACalculator : public PFCPositionCalculatorBase {
0022 public:
0023   Cluster3DPCACalculator(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0024       : PFCPositionCalculatorBase(conf, cc),
0025         updateTiming_(conf.getParameter<bool>("updateTiming")),
0026         pca_(new TPrincipal(3, "D")) {}
0027   Cluster3DPCACalculator(const Cluster3DPCACalculator&) = delete;
0028   Cluster3DPCACalculator& operator=(const Cluster3DPCACalculator&) = delete;
0029
0030   void calculateAndSetPosition(reco::PFCluster&, const HcalPFCuts*) override;
0031   void calculateAndSetPositions(reco::PFClusterCollection&, const HcalPFCuts*) override;
0032
0033 private:
0034   const bool updateTiming_;
0035   std::unique_ptr<TPrincipal> pca_;
0036
0037   void showerParameters(const reco::PFCluster&, math::XYZPoint&, math::XYZVector&);
0038
0039   void calculateAndSetPositionActual(reco::PFCluster&);
0040 };
0041
0042 DEFINE_EDM_PLUGIN(PFCPositionCalculatorFactory, Cluster3DPCACalculator, "Cluster3DPCACalculator");
0043
0044 void Cluster3DPCACalculator::calculateAndSetPosition(reco::PFCluster& cluster, const HcalPFCuts* cuts) {
0045   pca_ = std::make_unique<TPrincipal>(3, "D");
0046   calculateAndSetPositionActual(cluster);
0047 }
0048
0049 void Cluster3DPCACalculator::calculateAndSetPositions(reco::PFClusterCollection& clusters, const HcalPFCuts* cuts) {
0050   for (reco::PFCluster& cluster : clusters) {
0051     pca_ = std::make_unique<TPrincipal>(3, "D");
0052     calculateAndSetPositionActual(cluster);
0053   }
0054 }
0055
0056 void Cluster3DPCACalculator::calculateAndSetPositionActual(reco::PFCluster& cluster) {
0057   if (!cluster.seed()) {
0058     throw cms::Exception("ClusterWithNoSeed") << " Found a cluster with no seed: " << cluster;
0059   }
0060   double cl_energy = 0;
0061   double max_e = 0.0;
0062   double avg_time = 0.0;
0063   double time_norm = 0.0;
0064   PFLayer::Layer max_e_layer = PFLayer::NONE;
0065   reco::PFRecHitRef refseed;
0066   double pcavars[3];
0067
0068   for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
0069     const reco::PFRecHitRef& refhit = rhf.recHitRef();
0070     double rh_energy = refhit->energy();
0071     double rh_time = refhit->time();
0072     cl_energy += rh_energy * rhf.fraction();
0073     if (rh_time > 0.0) {  // time == -1 means no measurement
0074       // all times are offset by one nanosecond in digitizer
0075       // remove that here so all times of flight
0076       // are with respect to (0,0,0)
0077       avg_time += (rh_time - 1.0);
0078       time_norm += 1.0;
0079     }
0080     if (rh_energy > max_e) {
0081       max_e = rh_energy;
0082       max_e_layer = rhf.recHitRef()->layer();
0083     }
0084     if (refhit->detId() == cluster.seed())
0085       refseed = refhit;
0086     const double rh_fraction = rhf.fraction();
0087     rh_energy = refhit->energy() * rh_fraction;
0088     if (edm::isNotFinite(rh_energy)) {
0089       //temporarily changed exception to warning
0090       //      throw cms::Exception("PFClusterAlgo")
0091       edm::LogWarning("PFClusterAlgo") << "rechit " << refhit->detId() << " has a NaN energy... "
0092                                        << "The input of the particle flow clustering seems to be corrupted.";
0093       continue;
0094     }
0095     pcavars[0] = refhit->position().x();
0096     pcavars[1] = refhit->position().y();
0097     pcavars[2] = refhit->position().z();
0098     int nhit = int(rh_energy * 100);  // put rec_hit energy in units of 10 MeV
0099
0100     for (int i = 0; i < nhit; ++i) {
0101       pca_->AddRow(pcavars);
0102     }
0103   }
0104   cluster.setEnergy(cl_energy);
0105   cluster.setLayer(max_e_layer);
0106   // calculate the position
0107
0108   pca_->MakePrincipals();
0109   const TVectorD& means = *(pca_->GetMeanValues());
0110   const TMatrixD& eigens = *(pca_->GetEigenVectors());
0111
0112   math::XYZPoint barycenter(means[0], means[1], means[2]);
0113   math::XYZVector axis(eigens(0, 0), eigens(1, 0), eigens(2, 0));
0114
0115   if (time_norm > 0.0) {
0116     avg_time = avg_time / time_norm;
0117   } else {
0118     avg_time = std::numeric_limits<double>::min();
0119   }
0120
0121   if (axis.z() * barycenter.z() < 0.0) {
0122     axis = math::XYZVector(-eigens(0, 0), -eigens(1, 0), -eigens(2, 0));
0123   }
0124
0125   if (updateTiming_) {
0126     cluster.setTime(avg_time);
0127   }
0128   cluster.setPosition(barycenter);
0129   cluster.calculatePositionREP();
0130 }