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 #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   // for quick indexing back to hit energy
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     // Skip empty tracksters
0075     if (tst.vertices().empty()) {
0076       continue;
0077     }
0078     // Filter using trackster PID
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         if (cell.first.det() == DetId::Ecal || cell.first.det() == DetId::Hcal) {
0106           continue;
0107         }
0108         hitsAndFractions.emplace_back(cell.first, cell.second * fraction);
0109       }
0110     });
0111
0112     for (const auto& hAndF : hitsAndFractions) {
0113       auto itr = detIdToIndex.find(hAndF.first);
0114       if (itr == detIdToIndex.end()) {
0115         continue;  // hit wasn't saved in reco
0116       }
0117       auto ref = makeRefhit(input, itr->second);
0118       assert(ref->detId() == hAndF.first.rawId());
0119       const double hit_energy = hAndF.second * ref->energy();
0120       energy += hit_energy;
0121       back.addRecHitFraction(reco::PFRecHitFraction(ref, hAndF.second));
0122       // TODO: the following logic to identify the seed of a cluster
0123       // could be appropriate for the Run2 Ecal Calorimetric
0124       // detector, but could be wrong for the HGCal one. This has to
0125       // be reviewd.
0126       if (hit_energy > highest_energy || highest_energy == 0.0) {
0127         highest_energy = hit_energy;
0128         seed = ref->detId();
0129       }
0130     }  // end of hitsAndFractions
0131
0132     if (!back.hitsAndFractions().empty()) {
0133       back.setSeed(seed);
0134       back.setEnergy(energy);
0135       back.setCorrectedEnergy(energy);
0136     } else {
0137       back.setSeed(-1);
0138       back.setEnergy(0.f);
0139     }
0140   }  // end of loop over hgcalTracksters (3D)
0141 }