EgammaMCTools/plugins/PFClusterMatchedToPhotonsSelector.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CommonTools/RecoAlgos
0004 // Class:      PFClusterMatchedToPhotonsSelector
0005 //
0006 /**\class PFClusterMatchedToPhotonsSelector PFClusterMatchedToPhotonsSelector.cc CommonTools/RecoAlgos/plugins/PFClusterMatchedToPhotonsSelector.cc
0007
0008  Description: Matches ECAL PF clusters to photons that do not convert
0009
0010  Implementation:
0011
0012 */
0013 //
0014 // Original Author:  RCLSA
0015 //         Created:  Wed, 22 Mar 2017 18:01:40 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <iostream>
0022
0023 // user include files
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_;  // genParticles
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   //now do what ever initialization is needed
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     // Preselect PFclusters
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 //define this as a plug-in
0235 DEFINE_FWK_MODULE(PFClusterMatchedToPhotonsSelector);