TICL/plugins/MultiClustersFromTrackstersProducer.cc

0001 // user include files
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/stream/EDProducer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0009
0010 #include "DataFormats/ParticleFlowReco/interface/HGCalMultiCluster.h"
0011 #include "DataFormats/HGCalReco/interface/Trackster.h"
0012
0013 #include <algorithm>
0014 #include <array>
0015 #include <string>
0016 #include <vector>
0017
0018 class MultiClustersFromTrackstersProducer : public edm::stream::EDProducer<> {
0019 public:
0020   MultiClustersFromTrackstersProducer(const edm::ParameterSet&);
0021   ~MultiClustersFromTrackstersProducer() override {}
0022   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0023
0024   void produce(edm::Event&, const edm::EventSetup&) override;
0025
0026 private:
0027   edm::EDGetTokenT<std::vector<reco::CaloCluster>> layer_clusters_token_;
0028   edm::EDGetTokenT<std::vector<ticl::Trackster>> tracksters_token_;
0029 };
0030
0031 DEFINE_FWK_MODULE(MultiClustersFromTrackstersProducer);
0032
0033 MultiClustersFromTrackstersProducer::MultiClustersFromTrackstersProducer(const edm::ParameterSet& ps)
0034     : layer_clusters_token_(consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("LayerClusters"))),
0035       tracksters_token_(consumes<std::vector<ticl::Trackster>>(ps.getParameter<edm::InputTag>("Tracksters"))) {
0036   produces<std::vector<reco::HGCalMultiCluster>>();
0037 }
0038
0039 void MultiClustersFromTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0040   // hgcalMultiClusters
0041   edm::ParameterSetDescription desc;
0042   desc.add<edm::InputTag>("Tracksters", edm::InputTag("Tracksters", "TrackstersByCA"));
0043   desc.add<edm::InputTag>("LayerClusters", edm::InputTag("hgcalMergeLayerClusters"));
0044   desc.addUntracked<unsigned int>("verbosity", 3);
0045   descriptions.add("multiClustersFromTrackstersProducer", desc);
0046 }
0047
0048 void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0049   auto multiclusters = std::make_unique<std::vector<reco::HGCalMultiCluster>>();
0050   edm::Handle<std::vector<ticl::Trackster>> tracksterHandle;
0051   evt.getByToken(tracksters_token_, tracksterHandle);
0052
0053   edm::Handle<std::vector<reco::CaloCluster>> layer_clustersHandle;
0054   evt.getByToken(layer_clusters_token_, layer_clustersHandle);
0055
0056   auto const& tracksters = *tracksterHandle;
0057   auto const& layerClusters = *layer_clustersHandle;
0058
0059   edm::PtrVector<reco::BasicCluster> clusterPtrs;
0060   for (unsigned i = 0; i < layerClusters.size(); ++i) {
0061     edm::Ptr<reco::BasicCluster> ptr(layer_clustersHandle, i);
0062     clusterPtrs.push_back(ptr);
0063   }
0064
0065   std::for_each(std::begin(tracksters), std::end(tracksters), [&](auto const& trackster) {
0066     // Create an empty multicluster if the trackster has no layer clusters.
0067     // This could happen when a seed leads to no trackster and a dummy one is produced.
0068
0069     std::array<double, 3> baricenter{{0., 0., 0.}};
0070     double total_weight = 0.;
0071     reco::HGCalMultiCluster temp;
0072     int counter = 0;
0073     if (!trackster.vertices().empty()) {
0074       std::for_each(std::begin(trackster.vertices()), std::end(trackster.vertices()), [&](unsigned int idx) {
0075         temp.push_back(clusterPtrs[idx]);
0076         auto fraction = 1.f / trackster.vertex_multiplicity(counter++);
0077         for (auto const& cell : clusterPtrs[idx]->hitsAndFractions()) {
0078           temp.addHitAndFraction(cell.first, cell.second * fraction);
0079         }
0080         auto weight = clusterPtrs[idx]->energy() * fraction;
0081         total_weight += weight;
0082         baricenter[0] += clusterPtrs[idx]->x() * weight;
0083         baricenter[1] += clusterPtrs[idx]->y() * weight;
0084         baricenter[2] += clusterPtrs[idx]->z() * weight;
0085       });
0086       std::transform(
0087           std::begin(baricenter), std::end(baricenter), std::begin(baricenter), [&total_weight](double val) -> double {
0088             return val / total_weight;
0089           });
0090     }
0091     temp.setEnergy(total_weight);
0092     temp.setCorrectedEnergy(total_weight);
0093     temp.setPosition(math::XYZPoint(baricenter[0], baricenter[1], baricenter[2]));
0094     temp.setAlgoId(reco::CaloCluster::hgcal_em);
0095     temp.setTime(trackster.time(), trackster.timeError());
0096     multiclusters->push_back(temp);
0097   });
0098
0099   evt.put(std::move(multiclusters));
0100 }