File indexing completed on 2021-11-26 23:20:36
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
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
0071 if (tst.vertices().empty()) {
0072 continue;
0073 }
0074
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;
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
0116
0117
0118
0119 if (hit_energy > highest_energy || highest_energy == 0.0) {
0120 highest_energy = hit_energy;
0121 seed = ref->detId();
0122 }
0123 }
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 }
0134 }