File indexing completed on 2023-08-06 22:43:31
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <iostream>
0022
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/stream/EDProducer.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include "DataFormats/Common/interface/ValueMap.h"
0031 #include "DataFormats/Common/interface/Ref.h"
0032 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0033 #include "DataFormats/Math/interface/deltaR.h"
0034 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0035 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0036 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0037 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0038 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0039 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0040 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0041
0042 typedef reco::PFCluster::EEtoPSAssociation::value_type EEPSPair;
0043 bool sortByKey(const EEPSPair& a, const EEPSPair& b) { return a.first < b.first; }
0044
0045 class PFClusterMatchedToPhotonsSelector : public edm::stream::EDProducer<> {
0046 public:
0047 PFClusterMatchedToPhotonsSelector(const edm::ParameterSet&);
0048 static void fillDescriptions(edm::ConfigurationDescriptions&);
0049
0050 private:
0051 void produce(edm::Event&, const edm::EventSetup&) override;
0052
0053 edm::EDGetTokenT<reco::GenParticleCollection> genParticleToken_;
0054 edm::EDGetTokenT<reco::PFClusterCollection> particleFlowClusterECALToken_;
0055 edm::EDGetTokenT<reco::PFCluster::EEtoPSAssociation> associationToken_;
0056 edm::EDGetTokenT<TrackingParticleCollection> trackingParticleToken_;
0057 edm::EDGetTokenT<EcalRecHitCollection> recHitsEB_;
0058 edm::EDGetTokenT<EcalRecHitCollection> recHitsEE_;
0059
0060 const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_;
0061
0062 double matchMaxDR2_;
0063 double matchMaxDEDR2_;
0064
0065 double volumeZ_EB_;
0066 double volumeRadius_EB_;
0067 double volumeZ_EE_;
0068 };
0069
0070 PFClusterMatchedToPhotonsSelector::PFClusterMatchedToPhotonsSelector(const edm::ParameterSet& iConfig)
0071 : ecalClusterToolsESGetTokens_{consumesCollector()} {
0072
0073 particleFlowClusterECALToken_ =
0074 consumes<reco::PFClusterCollection>(iConfig.getParameter<edm::InputTag>("pfClustersTag"));
0075 associationToken_ =
0076 consumes<reco::PFCluster::EEtoPSAssociation>(iConfig.getParameter<edm::InputTag>("pfClustersTag"));
0077 trackingParticleToken_ =
0078 consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticleTag"));
0079 genParticleToken_ = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticleTag"));
0080 recHitsEB_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsEBLabel"));
0081 recHitsEE_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsEELabel"));
0082
0083 matchMaxDR2_ = iConfig.getParameter<double>("maxDR2");
0084 matchMaxDEDR2_ = iConfig.getParameter<double>("maxDEDR2");
0085 volumeZ_EB_ = iConfig.getParameter<double>("volumeZ_EB");
0086 volumeRadius_EB_ = iConfig.getParameter<double>("volumeRadius_EB");
0087 volumeZ_EE_ = iConfig.getParameter<double>("volumeZ_EE");
0088
0089 produces<reco::PFClusterCollection>();
0090 produces<edm::ValueMap<reco::GenParticleRef> >();
0091 produces<edm::ValueMap<int> >();
0092 produces<edm::ValueMap<float> >("PS1");
0093 produces<edm::ValueMap<float> >("PS2");
0094 }
0095
0096 void PFClusterMatchedToPhotonsSelector::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0097 edm::ParameterSetDescription desc;
0098 desc.add<edm::InputTag>("pfClustersTag", edm::InputTag("particleFlowClusterECAL"));
0099 desc.add<edm::InputTag>("trackingParticleTag", edm::InputTag("mix", "MergedTrackTruth"));
0100 desc.add<edm::InputTag>("genParticleTag", edm::InputTag("genParticles"));
0101 desc.add<edm::InputTag>("recHitsEBLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0102 desc.add<edm::InputTag>("recHitsEELabel", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0103 desc.add<double>("maxDR2", 0.1 * 0.1);
0104 desc.add<double>("maxDEDR2", 0.5 * 0.5);
0105 desc.add<double>("volumeZ_EB", 304.5);
0106 desc.add<double>("volumeRadius_EB", 123.8);
0107 desc.add<double>("volumeZ_EE", 317.0);
0108 descriptions.add("pfClusterMatchedToPhotonsSelector", desc);
0109 }
0110
0111 void PFClusterMatchedToPhotonsSelector::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0112 edm::Handle<reco::PFClusterCollection> particleFlowClusterECALHandle_;
0113 edm::Handle<reco::PFCluster::EEtoPSAssociation> associationHandle_;
0114 edm::Handle<TrackingParticleCollection> trackingParticleHandle_;
0115 edm::Handle<reco::GenParticleCollection> genParticleHandle_;
0116
0117 iEvent.getByToken(particleFlowClusterECALToken_, particleFlowClusterECALHandle_);
0118 iEvent.getByToken(trackingParticleToken_, trackingParticleHandle_);
0119 iEvent.getByToken(genParticleToken_, genParticleHandle_);
0120 iEvent.getByToken(associationToken_, associationHandle_);
0121
0122 std::unique_ptr<reco::PFClusterCollection> out = std::make_unique<reco::PFClusterCollection>();
0123 std::unique_ptr<edm::ValueMap<reco::GenParticleRef> > genmatching_out =
0124 std::make_unique<edm::ValueMap<reco::GenParticleRef> >();
0125 std::unique_ptr<edm::ValueMap<int> > clustersize_out = std::make_unique<edm::ValueMap<int> >();
0126 std::unique_ptr<edm::ValueMap<float> > energyPS1_out = std::make_unique<edm::ValueMap<float> >();
0127 std::unique_ptr<edm::ValueMap<float> > energyPS2_out = std::make_unique<edm::ValueMap<float> >();
0128
0129 std::vector<reco::GenParticleRef> genmatching;
0130 std::vector<int> clustersize;
0131 std::vector<float> energyPS1;
0132 std::vector<float> energyPS2;
0133
0134 EcalClusterLazyTools lazyTool(iEvent, ecalClusterToolsESGetTokens_.get(iSetup), recHitsEB_, recHitsEE_);
0135
0136 for (size_t iP = 0; iP < particleFlowClusterECALHandle_->size(); iP++) {
0137 auto&& pfCluster = particleFlowClusterECALHandle_->at(iP);
0138 bool isMatched = false;
0139 reco::GenParticleRef::key_type matchedKey;
0140
0141
0142 if (pfCluster.energy() <= 0) {
0143 continue;
0144 }
0145
0146 for (auto&& trackingParticle : *trackingParticleHandle_) {
0147 if (trackingParticle.pdgId() != 22)
0148 continue;
0149 if (trackingParticle.status() != 1)
0150 continue;
0151 matchedKey = trackingParticle.genParticles().at(0).key();
0152
0153 float dR2 = reco::deltaR2(trackingParticle, pfCluster.position());
0154 if (dR2 > matchMaxDR2_)
0155 continue;
0156 float dE = 1. - trackingParticle.genParticles().at(0)->energy() / pfCluster.energy();
0157 if ((dR2 + dE * dE) > matchMaxDEDR2_)
0158 continue;
0159
0160 bool isConversion = false;
0161 for (auto&& vertRef : trackingParticle.decayVertices()) {
0162 if (vertRef->position().rho() > volumeRadius_EB_ && std::abs(vertRef->position().z()) < volumeZ_EB_)
0163 continue;
0164 if (std::abs(vertRef->position().z()) > volumeZ_EE_)
0165 continue;
0166
0167 for (auto&& tpRef : vertRef->daughterTracks()) {
0168 if (std::abs(tpRef->pdgId()) == 11)
0169 isConversion = true;
0170 break;
0171 }
0172 if (isConversion)
0173 break;
0174 }
0175 if (isConversion)
0176 continue;
0177
0178 isMatched = true;
0179 break;
0180 }
0181
0182 if (isMatched) {
0183 out->push_back(pfCluster);
0184 double ePS1 = 0, ePS2 = 0;
0185 if (!(pfCluster.layer() == PFLayer::ECAL_BARREL)) {
0186 auto ee_key_val = std::make_pair(iP, edm::Ptr<reco::PFCluster>());
0187 const auto clustops = std::equal_range(
0188 associationHandle_.product()->begin(), associationHandle_.product()->end(), ee_key_val, sortByKey);
0189 for (auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) {
0190 edm::Ptr<reco::PFCluster> psclus(i_ps->second);
0191 switch (psclus->layer()) {
0192 case PFLayer::PS1:
0193 ePS1 += psclus->energy();
0194 break;
0195 case PFLayer::PS2:
0196 ePS2 += psclus->energy();
0197 break;
0198 default:
0199 break;
0200 }
0201 }
0202 }
0203
0204 genmatching.push_back(edm::Ref<reco::GenParticleCollection>(genParticleHandle_, matchedKey));
0205 clustersize.push_back(lazyTool.n5x5(pfCluster));
0206 energyPS1.push_back(ePS1);
0207 energyPS2.push_back(ePS2);
0208 }
0209 }
0210
0211 edm::OrphanHandle<reco::PFClusterCollection> pfClusterHandle = iEvent.put(std::move(out));
0212
0213 edm::ValueMap<reco::GenParticleRef>::Filler mapFiller(*genmatching_out);
0214 mapFiller.insert(pfClusterHandle, genmatching.begin(), genmatching.end());
0215 mapFiller.fill();
0216 iEvent.put(std::move(genmatching_out));
0217
0218 edm::ValueMap<int>::Filler mapFiller_int(*clustersize_out);
0219 mapFiller_int.insert(pfClusterHandle, clustersize.begin(), clustersize.end());
0220 mapFiller_int.fill();
0221 iEvent.put(std::move(clustersize_out));
0222
0223 edm::ValueMap<float>::Filler mapFiller_energyPS1(*energyPS1_out);
0224 mapFiller_energyPS1.insert(pfClusterHandle, energyPS1.begin(), energyPS1.end());
0225 mapFiller_energyPS1.fill();
0226 iEvent.put(std::move(energyPS1_out), "PS1");
0227
0228 edm::ValueMap<float>::Filler mapFiller_energyPS2(*energyPS2_out);
0229 mapFiller_energyPS2.insert(pfClusterHandle, energyPS2.begin(), energyPS2.end());
0230 mapFiller_energyPS2.fill();
0231 iEvent.put(std::move(energyPS2_out), "PS2");
0232 }
0233
0234
0235 DEFINE_FWK_MODULE(PFClusterMatchedToPhotonsSelector);