File indexing completed on 2023-10-25 09:59:08
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
0060 auto gsfPfRecTracks = edm::makeValid(event.getHandle(gsfPfRecTracks_));
0061
0062
0063 reco::PFClusterCollection const& ecalClusters = event.get(ecalClusters_);
0064
0065
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
0071 std::vector<reco::SuperClusterRef> superClustersValueMap;
0072
0073
0074 auto caloClusters = std::make_unique<reco::CaloClusterCollection>(reco::CaloClusterCollection());
0075 caloClusters->reserve(ecalClusters.size());
0076
0077
0078 std::vector<std::vector<int> > cluster_idx;
0079 cluster_idx.reserve(gsfPfRecTracks->size());
0080
0081
0082 std::vector<std::vector<int> > pfcluster_idx;
0083 pfcluster_idx.reserve(gsfPfRecTracks->size());
0084
0085
0086 std::vector<std::vector<float> > cluster_dr2min;
0087 cluster_dr2min.reserve(gsfPfRecTracks->size());
0088
0089
0090 std::vector<std::vector<const reco::PFTrajectoryPoint*> > points;
0091 points.reserve(gsfPfRecTracks->size());
0092 for (auto const& trk : *gsfPfRecTracks) {
0093
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
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
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
0110 for (size_t iclu = 0; iclu < ecalClusters.size(); ++iclu) {
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) {
0114 for (size_t jpoint = 0; jpoint < points[ipoint].size(); ++jpoint) {
0115 if (points[ipoint][jpoint]->isValid()) {
0116 float dr2 = reco::deltaR2(ecalClusters[iclu], points[ipoint][jpoint]->positionREP());
0117 if (dr2 < dr2min) {
0118
0119 dr2min = dr2;
0120 point = std::make_pair(ipoint, jpoint);
0121 }
0122 }
0123 }
0124 }
0125 if (point.first >= 0 && point.second >= 0 &&
0126 dr2min < cluster_dr2min[point.first][point.second]) {
0127
0128 caloClusters->push_back(ecalClusters[iclu]);
0129
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
0137 const edm::OrphanHandle<reco::CaloClusterCollection>& caloClustersH = event.put(std::move(caloClusters));
0138
0139
0140 for (size_t itrk = 0; itrk < gsfPfRecTracks->size(); ++itrk) {
0141
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
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
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_ ||
0171 index == cluster_idx[itrk][ipoint])) {
0172 continue;
0173 }
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
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();
0201 superClusters->push_back(sc);
0202
0203
0204 superClustersValueMap.push_back(reco::SuperClusterRef(superClustersRefProd, superClusters->size() - 1));
0205 } else {
0206 superClustersValueMap.push_back(reco::SuperClusterRef(superClustersRefProd.id()));
0207 }
0208
0209 }
0210
0211
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);