File indexing completed on 2024-04-06 12:27:24
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 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0008 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0009
0010 class PFClusterFromHGCalTrackster : public InitialClusteringStepBase {
0011 public:
0012 PFClusterFromHGCalTrackster(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0013 : InitialClusteringStepBase(conf, cc) {
0014 filterByTracksterPID_ = conf.getParameter<bool>("filterByTracksterPID");
0015 filterByTracksterIteration_ = conf.getParameter<bool>("filterByTracksterIteration");
0016 pid_threshold_ = conf.getParameter<double>("pid_threshold");
0017 filter_on_categories_ = conf.getParameter<std::vector<int> >("filter_on_categories");
0018 filter_on_iterations_ = conf.getParameter<std::vector<int> >("filter_on_iterations");
0019
0020 tracksterToken_ = cc.consumes<std::vector<ticl::Trackster> >(conf.getParameter<edm::InputTag>("tracksterSrc"));
0021 clusterToken_ = cc.consumes<reco::CaloClusterCollection>(conf.getParameter<edm::InputTag>("clusterSrc"));
0022 }
0023
0024 ~PFClusterFromHGCalTrackster() override {}
0025 PFClusterFromHGCalTrackster(const PFClusterFromHGCalTrackster&) = delete;
0026 PFClusterFromHGCalTrackster& operator=(const PFClusterFromHGCalTrackster&) = delete;
0027
0028 void updateEvent(const edm::Event&) final;
0029
0030 void buildClusters(const edm::Handle<reco::PFRecHitCollection>&,
0031 const std::vector<bool>&,
0032 const std::vector<bool>&,
0033 reco::PFClusterCollection&,
0034 const HcalPFCuts*) override;
0035
0036 private:
0037 bool filterByTracksterPID_;
0038 bool filterByTracksterIteration_;
0039 float pid_threshold_;
0040 std::vector<int> filter_on_categories_;
0041 std::vector<int> filter_on_iterations_;
0042
0043 edm::EDGetTokenT<std::vector<ticl::Trackster> > tracksterToken_;
0044 edm::Handle<std::vector<ticl::Trackster> > trackstersH_;
0045
0046 edm::EDGetTokenT<reco::CaloClusterCollection> clusterToken_;
0047 edm::Handle<reco::CaloClusterCollection> clusterH_;
0048 };
0049
0050 DEFINE_EDM_PLUGIN(InitialClusteringStepFactory, PFClusterFromHGCalTrackster, "PFClusterFromHGCalTrackster");
0051
0052 void PFClusterFromHGCalTrackster::updateEvent(const edm::Event& ev) {
0053 ev.getByToken(tracksterToken_, trackstersH_);
0054 ev.getByToken(clusterToken_, clusterH_);
0055 }
0056
0057 void PFClusterFromHGCalTrackster::buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
0058 const std::vector<bool>& rechitMask,
0059 const std::vector<bool>& seedable,
0060 reco::PFClusterCollection& output,
0061 const HcalPFCuts*) {
0062 auto const& hits = *input;
0063
0064 const auto& tracksters = *trackstersH_;
0065 const auto& clusters = *clusterH_;
0066
0067
0068 std::unordered_map<uint32_t, size_t> detIdToIndex(hits.size());
0069 for (uint32_t i = 0; i < hits.size(); ++i) {
0070 detIdToIndex[hits[i].detId()] = i;
0071 }
0072
0073 for (const auto& tst : tracksters) {
0074
0075 if (tst.vertices().empty()) {
0076 continue;
0077 }
0078
0079 if (filterByTracksterPID_) {
0080 float probTotal = 0.0f;
0081 for (int cat : filter_on_categories_) {
0082 probTotal += tst.id_probabilities(cat);
0083 }
0084 if (probTotal < pid_threshold_) {
0085 continue;
0086 }
0087 }
0088
0089 if (filterByTracksterIteration_ &&
0090 std::find(filter_on_iterations_.begin(), filter_on_iterations_.end(), tst.ticlIteration()) ==
0091 filter_on_iterations_.end()) {
0092 continue;
0093 }
0094
0095 DetId seed;
0096 double energy = 0.0, highest_energy = 0.0;
0097 output.emplace_back();
0098 reco::PFCluster& back = output.back();
0099
0100 std::vector<std::pair<DetId, float> > hitsAndFractions;
0101 int iLC = 0;
0102 std::for_each(std::begin(tst.vertices()), std::end(tst.vertices()), [&](unsigned int lcId) {
0103 const auto fraction = 1.f / tst.vertex_multiplicity(iLC++);
0104 for (const auto& cell : clusters[lcId].hitsAndFractions()) {
0105 hitsAndFractions.emplace_back(cell.first, cell.second * fraction);
0106 }
0107 });
0108
0109 for (const auto& hAndF : hitsAndFractions) {
0110 auto itr = detIdToIndex.find(hAndF.first);
0111 if (itr == detIdToIndex.end()) {
0112 continue;
0113 }
0114 auto ref = makeRefhit(input, itr->second);
0115 assert(ref->detId() == hAndF.first.rawId());
0116 const double hit_energy = hAndF.second * ref->energy();
0117 energy += hit_energy;
0118 back.addRecHitFraction(reco::PFRecHitFraction(ref, hAndF.second));
0119
0120
0121
0122
0123 if (hit_energy > highest_energy || highest_energy == 0.0) {
0124 highest_energy = hit_energy;
0125 seed = ref->detId();
0126 }
0127 }
0128
0129 if (!back.hitsAndFractions().empty()) {
0130 back.setSeed(seed);
0131 back.setEnergy(energy);
0132 back.setCorrectedEnergy(energy);
0133 } else {
0134 back.setSeed(0);
0135 back.setEnergy(0.f);
0136 }
0137 }
0138 }