EgammaElectronProducers/plugins/LowPtGsfElectronSCProducer.cc

0001 #include "CommonTools/ParticleFlow/interface/PFClusterWidthAlgo.h"
0002 #include "DataFormats/Common/interface/RefToPtr.h"
0003 #include "DataFormats/Common/interface/ValueMap.h"
0004 #include "DataFormats/Common/interface/ValidHandle.h"
0005 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0006 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0007 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0008 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0009 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0010 #include "DataFormats/Math/interface/deltaR.h"
0011 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
0012 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
0017 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrackFwd.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0020 #include "DataFormats/ParticleFlowReco/interface/PFTrajectoryPoint.h"
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/stream/EDProducer.h"
0023 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0024
0025 #include <iostream>
0026
0027 class LowPtGsfElectronSCProducer : public edm::stream::EDProducer<> {
0028 public:
0029   explicit LowPtGsfElectronSCProducer(const edm::ParameterSet&);
0030
0031   void produce(edm::Event&, const edm::EventSetup&) override;
0032
0033   static void fillDescriptions(edm::ConfigurationDescriptions&);
0034
0035 private:
0036   reco::PFClusterRef closestCluster(const reco::PFTrajectoryPoint& point,
0037                                     const edm::Handle<reco::PFClusterCollection>& clusters,
0038                                     std::vector<int>& matched);
0039
0040   const edm::EDGetTokenT<reco::GsfPFRecTrackCollection> gsfPfRecTracks_;
0041   const edm::EDGetTokenT<reco::PFClusterCollection> ecalClusters_;
0042   const double dr2_;
0043 };
0044
0045 ////////////////////////////////////////////////////////////////////////////////
0046 //
0047 LowPtGsfElectronSCProducer::LowPtGsfElectronSCProducer(const edm::ParameterSet& cfg)
0048     : gsfPfRecTracks_{consumes<reco::GsfPFRecTrackCollection>(cfg.getParameter<edm::InputTag>("gsfPfRecTracks"))},
0049       ecalClusters_{consumes<reco::PFClusterCollection>(cfg.getParameter<edm::InputTag>("ecalClusters"))},
0050       dr2_{cfg.getParameter<double>("MaxDeltaR2")} {
0051   produces<reco::CaloClusterCollection>();
0052   produces<reco::SuperClusterCollection>();
0053   produces<edm::ValueMap<reco::SuperClusterRef> >();
0054 }
0055
0056 ////////////////////////////////////////////////////////////////////////////////
0057 //
0058 void LowPtGsfElectronSCProducer::produce(edm::Event& event, const edm::EventSetup&) {
0059   // Input GsfPFRecTracks collection
0060   auto gsfPfRecTracks = edm::makeValid(event.getHandle(gsfPfRecTracks_));
0061
0062   // Input EcalClusters collection
0063   reco::PFClusterCollection const& ecalClusters = event.get(ecalClusters_);
0064
0065   // Output SuperClusters collection and getRefBeforePut
0066   auto superClusters = std::make_unique<reco::SuperClusterCollection>(reco::SuperClusterCollection());
0067   superClusters->reserve(gsfPfRecTracks->size());
0068   const reco::SuperClusterRefProd superClustersRefProd = event.getRefBeforePut<reco::SuperClusterCollection>();
0069
0070   // Output ValueMap container of GsfPFRecTrackRef index to SuperClusterRef
0071   std::vector<reco::SuperClusterRef> superClustersValueMap;
0072
0073   // Output CaloClusters collection
0074   auto caloClusters = std::make_unique<reco::CaloClusterCollection>(reco::CaloClusterCollection());
0075   caloClusters->reserve(ecalClusters.size());
0076
0077   // Index[GSF track][trajectory point] for "best" CaloCluster
0078   std::vector<std::vector<int> > cluster_idx;
0079   cluster_idx.reserve(gsfPfRecTracks->size());
0080
0081   // Index[GSF track][trajectory point] for "best" PFCluster
0082   std::vector<std::vector<int> > pfcluster_idx;
0083   pfcluster_idx.reserve(gsfPfRecTracks->size());
0084
0085   // dr2min[GSF track][trajectory point] for "best" CaloCluster
0086   std::vector<std::vector<float> > cluster_dr2min;
0087   cluster_dr2min.reserve(gsfPfRecTracks->size());
0088
0089   // Construct list of trajectory points from the GSF track and electron brems
0090   std::vector<std::vector<const reco::PFTrajectoryPoint*> > points;
0091   points.reserve(gsfPfRecTracks->size());
0092   for (auto const& trk : *gsfPfRecTracks) {
0093     // Extrapolated track
0094     std::vector<const reco::PFTrajectoryPoint*> traj;
0095     traj.reserve(trk.PFRecBrem().size() + 1);
0096     traj.push_back(&trk.extrapolatedPoint(reco::PFTrajectoryPoint::LayerType::ECALShowerMax));
0097     // Extrapolated brem trajectories
0098     for (auto const& brem : trk.PFRecBrem()) {
0099       traj.push_back(&brem.extrapolatedPoint(reco::PFTrajectoryPoint::LayerType::ECALShowerMax));
0100     }
0101     auto size = traj.size();
0102     points.push_back(std::move(traj));
0103     // Size containers
0104     cluster_idx.emplace_back(size, -1);
0105     pfcluster_idx.emplace_back(size, -1);
0106     cluster_dr2min.emplace_back(size, 1.e6);
0107   }
0108
0109   // For each cluster, find closest trajectory point ("global" arbitration at event level)
0110   for (size_t iclu = 0; iclu < ecalClusters.size(); ++iclu) {  // Cluster loop
0111     std::pair<int, int> point = std::make_pair(-1, -1);
0112     float dr2min = 1.e6;
0113     for (size_t ipoint = 0; ipoint < points.size(); ++ipoint) {            // GSF track loop
0114       for (size_t jpoint = 0; jpoint < points[ipoint].size(); ++jpoint) {  // Traj point loop
0115         if (points[ipoint][jpoint]->isValid()) {
0116           float dr2 = reco::deltaR2(ecalClusters[iclu], points[ipoint][jpoint]->positionREP());
0117           if (dr2 < dr2min) {
0118             // Store nearest point to this cluster
0119             dr2min = dr2;
0120             point = std::make_pair(ipoint, jpoint);
0121           }
0122         }
0123       }
0124     }
0125     if (point.first >= 0 && point.second >= 0 &&               // if this cluster is matched to a point ...
0126         dr2min < cluster_dr2min[point.first][point.second]) {  // ... and cluster is closest to the same point
0127       // Copy CaloCluster to new collection
0128       caloClusters->push_back(ecalClusters[iclu]);
0129       // Store cluster index for creation of SC later
0130       cluster_idx[point.first][point.second] = caloClusters->size() - 1;
0131       pfcluster_idx[point.first][point.second] = iclu;
0132       cluster_dr2min[point.first][point.second] = dr2min;
0133     }
0134   }
0135
0136   // Put CaloClusters in event and get orphan handle
0137   const edm::OrphanHandle<reco::CaloClusterCollection>& caloClustersH = event.put(std::move(caloClusters));
0138
0139   // Loop through GSF tracks
0140   for (size_t itrk = 0; itrk < gsfPfRecTracks->size(); ++itrk) {
0141     // Used to create SC
0142     float energy = 0.;
0143     float X = 0., Y = 0., Z = 0.;
0144     reco::CaloClusterPtr seed;
0145     reco::CaloClusterPtrVector clusters;
0146     std::vector<const reco::PFCluster*> barePtrs;
0147
0148     // Find closest match in dr2 from points associated to given track
0149     int index = -1;
0150     float dr2 = 1.e6;
0151     for (size_t ipoint = 0; ipoint < cluster_idx[itrk].size(); ++ipoint) {
0152       if (cluster_idx[itrk][ipoint] < 0) {
0153         continue;
0154       }
0155       if (cluster_dr2min[itrk][ipoint] < dr2) {
0156         dr2 = cluster_dr2min[itrk][ipoint];
0157         index = cluster_idx[itrk][ipoint];
0158       }
0159     }
0160
0161     // For each track, loop through points and use associated cluster
0162     for (size_t ipoint = 0; ipoint < cluster_idx[itrk].size(); ++ipoint) {
0163       if (cluster_idx[itrk][ipoint] < 0) {
0164         continue;
0165       }
0166       reco::CaloClusterPtr clu(caloClustersH, cluster_idx[itrk][ipoint]);
0167       if (clu.isNull()) {
0168         continue;
0169       }
0170       if (!(cluster_dr2min[itrk][ipoint] < dr2_ ||  // Require cluster to be closer than dr2_ ...
0171             index == cluster_idx[itrk][ipoint])) {
0172         continue;
0173       }  // ... unless it is the closest one ...
0174       if (seed.isNull()) {
0175         seed = clu;
0176       }
0177       clusters.push_back(clu);
0178       energy += clu->correctedEnergy();
0179       X += clu->position().X() * clu->correctedEnergy();
0180       Y += clu->position().Y() * clu->correctedEnergy();
0181       Z += clu->position().Z() * clu->correctedEnergy();
0182       auto index = pfcluster_idx[itrk][ipoint];
0183       if (index < static_cast<decltype(index)>(ecalClusters.size())) {
0184         barePtrs.push_back(&(ecalClusters[index]));
0185       }
0186     }
0187     X /= energy;
0188     Y /= energy;
0189     Z /= energy;
0190
0191     // Create SC
0192     if (seed.isNonnull()) {
0193       reco::SuperCluster sc(energy, math::XYZPoint(X, Y, Z));
0194       sc.setCorrectedEnergy(energy);
0195       sc.setSeed(seed);
0196       sc.setClusters(clusters);
0197       PFClusterWidthAlgo pfwidth(barePtrs);
0198       sc.setEtaWidth(pfwidth.pflowEtaWidth());
0199       sc.setPhiWidth(pfwidth.pflowPhiWidth());
0200       sc.rawEnergy();  // Cache the value of raw energy
0201       superClusters->push_back(sc);
0202
0203       // Populate ValueMap container
0204       superClustersValueMap.push_back(reco::SuperClusterRef(superClustersRefProd, superClusters->size() - 1));
0205     } else {
0206       superClustersValueMap.push_back(reco::SuperClusterRef(superClustersRefProd.id()));
0207     }
0208
0209   }  // GSF tracks
0210
0211   // Put SuperClusters in event
0212   event.put(std::move(superClusters));
0213
0214   auto ptr = std::make_unique<edm::ValueMap<reco::SuperClusterRef> >(edm::ValueMap<reco::SuperClusterRef>());
0215   edm::ValueMap<reco::SuperClusterRef>::Filler filler(*ptr);
0216   filler.insert(gsfPfRecTracks, superClustersValueMap.begin(), superClustersValueMap.end());
0217   filler.fill();
0218   event.put(std::move(ptr));
0219 }
0220
0221 ////////////////////////////////////////////////////////////////////////////////
0222 //
0223 reco::PFClusterRef LowPtGsfElectronSCProducer::closestCluster(const reco::PFTrajectoryPoint& point,
0224                                                               const edm::Handle<reco::PFClusterCollection>& clusters,
0225                                                               std::vector<int>& matched) {
0226   reco::PFClusterRef closest;
0227   if (point.isValid()) {
0228     float dr2min = dr2_;
0229     for (size_t ii = 0; ii < clusters->size(); ++ii) {
0230       if (std::find(matched.begin(), matched.end(), ii) == matched.end()) {
0231         float dr2 = reco::deltaR2(clusters->at(ii), point.positionREP());
0232         if (dr2 < dr2min) {
0233           closest = reco::PFClusterRef(clusters, ii);
0234           dr2min = dr2;
0235         }
0236       }
0237     }
0238     if (dr2min < (dr2_ - 1.e-6)) {
0239       matched.push_back(closest.index());
0240     }
0241   }
0242   return closest;
0243 }
0244
0245 //////////////////////////////////////////////////////////////////////////////////////////
0246 //
0247 void LowPtGsfElectronSCProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0248   edm::ParameterSetDescription desc;
0249   desc.add<edm::InputTag>("gsfPfRecTracks", edm::InputTag("lowPtGsfElePfGsfTracks"));
0250   desc.add<edm::InputTag>("ecalClusters", edm::InputTag("particleFlowClusterECAL"));
0251   desc.add<edm::InputTag>("hcalClusters", edm::InputTag("particleFlowClusterHCAL"));
0252   desc.add<double>("MaxDeltaR2", 0.5);
0253   descriptions.add("lowPtGsfElectronSuperClusters", desc);
0254 }
0255
0256 //////////////////////////////////////////////////////////////////////////////////////////
0257 //
0258 #include "FWCore/Framework/interface/MakerMacros.h"
0259 DEFINE_FWK_MODULE(LowPtGsfElectronSCProducer);