PFClusterProducer/plugins/PFClusterFromHGCalTrackster.cc

0001 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0002 #include "DataFormats/HGCalReco/interface/Trackster.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
0007
0008 class PFClusterFromHGCalTrackster : public InitialClusteringStepBase {
0009 public:
0010   PFClusterFromHGCalTrackster(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0011       : InitialClusteringStepBase(conf, cc) {
0012     filterByTracksterPID_ = conf.getParameter<bool>("filterByTracksterPID");
0013     filterByTracksterIteration_ = conf.getParameter<bool>("filterByTracksterIteration");
0014     pid_threshold_ = conf.getParameter<double>("pid_threshold");
0015     filter_on_categories_ = conf.getParameter<std::vector<int> >("filter_on_categories");
0016     filter_on_iterations_ = conf.getParameter<std::vector<int> >("filter_on_iterations");
0017
0018     tracksterToken_ = cc.consumes<std::vector<ticl::Trackster> >(conf.getParameter<edm::InputTag>("tracksterSrc"));
0019     clusterToken_ = cc.consumes<reco::CaloClusterCollection>(conf.getParameter<edm::InputTag>("clusterSrc"));
0020   }
0021
0022   ~PFClusterFromHGCalTrackster() override {}
0023   PFClusterFromHGCalTrackster(const PFClusterFromHGCalTrackster&) = delete;
0024   PFClusterFromHGCalTrackster& operator=(const PFClusterFromHGCalTrackster&) = delete;
0025
0026   void updateEvent(const edm::Event&) final;
0027
0028   void buildClusters(const edm::Handle<reco::PFRecHitCollection>&,
0029                      const std::vector<bool>&,
0030                      const std::vector<bool>&,
0031                      reco::PFClusterCollection&) override;
0032
0033 private:
0034   bool filterByTracksterPID_;
0035   bool filterByTracksterIteration_;
0036   float pid_threshold_;
0037   std::vector<int> filter_on_categories_;
0038   std::vector<int> filter_on_iterations_;
0039
0040   edm::EDGetTokenT<std::vector<ticl::Trackster> > tracksterToken_;
0041   edm::Handle<std::vector<ticl::Trackster> > trackstersH_;
0042
0043   edm::EDGetTokenT<reco::CaloClusterCollection> clusterToken_;
0044   edm::Handle<reco::CaloClusterCollection> clusterH_;
0045 };
0046
0047 DEFINE_EDM_PLUGIN(InitialClusteringStepFactory, PFClusterFromHGCalTrackster, "PFClusterFromHGCalTrackster");
0048
0049 void PFClusterFromHGCalTrackster::updateEvent(const edm::Event& ev) {
0050   ev.getByToken(tracksterToken_, trackstersH_);
0051   ev.getByToken(clusterToken_, clusterH_);
0052 }
0053
0054 void PFClusterFromHGCalTrackster::buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
0055                                                 const std::vector<bool>& rechitMask,
0056                                                 const std::vector<bool>& seedable,
0057                                                 reco::PFClusterCollection& output) {
0058   auto const& hits = *input;
0059
0060   const auto& tracksters = *trackstersH_;
0061   const auto& clusters = *clusterH_;
0062
0063   // for quick indexing back to hit energy
0064   std::unordered_map<uint32_t, size_t> detIdToIndex(hits.size());
0065   for (uint32_t i = 0; i < hits.size(); ++i) {
0066     detIdToIndex[hits[i].detId()] = i;
0067   }
0068
0069   for (const auto& tst : tracksters) {
0070     // Skip empty tracksters
0071     if (tst.vertices().empty()) {
0072       continue;
0073     }
0074     // Filter using trackster PID
0075     if (filterByTracksterPID_) {
0076       float probTotal = 0.0f;
0077       for (int cat : filter_on_categories_) {
0078         probTotal += tst.id_probabilities(cat);
0079       }
0080       if (probTotal < pid_threshold_) {
0081         continue;
0082       }
0083     }
0084
0085     if (filterByTracksterIteration_ &&
0086         std::find(filter_on_iterations_.begin(), filter_on_iterations_.end(), tst.ticlIteration()) ==
0087             filter_on_iterations_.end()) {
0088       continue;
0089     }
0090
0091     DetId seed;
0092     double energy = 0.0, highest_energy = 0.0;
0093     output.emplace_back();
0094     reco::PFCluster& back = output.back();
0095
0096     std::vector<std::pair<DetId, float> > hitsAndFractions;
0097     int iLC = 0;
0098     std::for_each(std::begin(tst.vertices()), std::end(tst.vertices()), [&](unsigned int lcId) {
0099       const auto fraction = 1.f / tst.vertex_multiplicity(iLC++);
0100       for (const auto& cell : clusters[lcId].hitsAndFractions()) {
0101         hitsAndFractions.emplace_back(cell.first, cell.second * fraction);
0102       }
0103     });
0104
0105     for (const auto& hAndF : hitsAndFractions) {
0106       auto itr = detIdToIndex.find(hAndF.first);
0107       if (itr == detIdToIndex.end()) {
0108         continue;  // hit wasn't saved in reco
0109       }
0110       auto ref = makeRefhit(input, itr->second);
0111       assert(ref->detId() == hAndF.first.rawId());
0112       const double hit_energy = hAndF.second * ref->energy();
0113       energy += hit_energy;
0114       back.addRecHitFraction(reco::PFRecHitFraction(ref, hAndF.second));
0115       // TODO: the following logic to identify the seed of a cluster
0116       // could be appropriate for the Run2 Ecal Calorimetric
0117       // detector, but could be wrong for the HGCal one. This has to
0118       // be reviewd.
0119       if (hit_energy > highest_energy || highest_energy == 0.0) {
0120         highest_energy = hit_energy;
0121         seed = ref->detId();
0122       }
0123     }  // end of hitsAndFractions
0124
0125     if (!back.hitsAndFractions().empty()) {
0126       back.setSeed(seed);
0127       back.setEnergy(energy);
0128       back.setCorrectedEnergy(energy);
0129     } else {
0130       back.setSeed(0);
0131       back.setEnergy(0.f);
0132     }
0133   }  // end of loop over hgcalTracksters (3D)
0134 }