EgammaIsolationAlgos/plugins/ParticleBasedIsoProducer.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "DataFormats/Common/interface/ValueMap.h"
0003 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0004 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0005 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0006 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0009 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0010 #include "FWCore/Framework/interface/stream/EDProducer.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "RecoEgamma/EgammaIsolationAlgos/interface/PFBlockBasedIsolation.h"
0015
0016 class ParticleBasedIsoProducer : public edm::stream::EDProducer<> {
0017 public:
0018   ParticleBasedIsoProducer(const edm::ParameterSet& conf);
0019
0020   void beginRun(edm::Run const& r, edm::EventSetup const& es) override;
0021   void produce(edm::Event& e, const edm::EventSetup& c) override;
0022
0023   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0024
0025 private:
0026   edm::ParameterSet conf_;
0027   std::string photonCollection_;
0028   std::string electronCollection_;
0029
0030   edm::InputTag photonProducer_;
0031   edm::InputTag photonTmpProducer_;
0032
0033   edm::InputTag electronProducer_;
0034   edm::InputTag electronTmpProducer_;
0035
0036   edm::EDGetTokenT<reco::PhotonCollection> photonProducerT_;
0037   edm::EDGetTokenT<reco::PhotonCollection> photonTmpProducerT_;
0038   edm::EDGetTokenT<reco::GsfElectronCollection> electronProducerT_;
0039   edm::EDGetTokenT<reco::GsfElectronCollection> electronTmpProducerT_;
0040   edm::EDGetTokenT<reco::PFCandidateCollection> pfEgammaCandidates_;
0041   edm::EDGetTokenT<reco::PFCandidateCollection> pfCandidates_;
0042   edm::EDGetTokenT<edm::ValueMap<reco::PhotonRef>> valMapPFCandToPhoton_;
0043   edm::EDGetTokenT<edm::ValueMap<reco::GsfElectronRef>> valMapPFCandToEle_;
0044
0045   std::string valueMapPFCandPhoton_;
0046   std::string valueMapPhoPFCandIso_;
0047   std::string valueMapPFCandEle_;
0048   std::string valueMapElePFCandIso_;
0049
0050   std::unique_ptr<PFBlockBasedIsolation> thePFBlockBasedIsolation_;
0051 };
0052
0053 #include "FWCore/Framework/interface/MakerMacros.h"
0054 DEFINE_FWK_MODULE(ParticleBasedIsoProducer);
0055
0056 ParticleBasedIsoProducer::ParticleBasedIsoProducer(const edm::ParameterSet& conf) : conf_(conf) {
0057   photonTmpProducer_ = conf_.getParameter<edm::InputTag>("photonTmpProducer");
0058   photonProducer_ = conf_.getParameter<edm::InputTag>("photonProducer");
0059   electronProducer_ = conf_.getParameter<edm::InputTag>("electronProducer");
0060   electronTmpProducer_ = conf_.getParameter<edm::InputTag>("electronTmpProducer");
0061
0062   photonProducerT_ = consumes<reco::PhotonCollection>(photonProducer_);
0063
0064   photonTmpProducerT_ = consumes<reco::PhotonCollection>(photonTmpProducer_);
0065
0066   electronProducerT_ = consumes<reco::GsfElectronCollection>(electronProducer_);
0067
0068   electronTmpProducerT_ = consumes<reco::GsfElectronCollection>(electronTmpProducer_);
0069
0070   pfCandidates_ = consumes<reco::PFCandidateCollection>(conf_.getParameter<edm::InputTag>("pfCandidates"));
0071
0072   pfEgammaCandidates_ = consumes<reco::PFCandidateCollection>(conf_.getParameter<edm::InputTag>("pfEgammaCandidates"));
0073
0074   valueMapPFCandPhoton_ = conf_.getParameter<std::string>("valueMapPhoToEG");
0075   valueMapPFCandEle_ = conf_.getParameter<std::string>("valueMapEleToEG");
0076
0077   valMapPFCandToPhoton_ = consumes<edm::ValueMap<reco::PhotonRef>>({"gedPhotonsTmp", valueMapPFCandPhoton_});
0078
0079   valMapPFCandToEle_ =
0080       consumes<edm::ValueMap<reco::GsfElectronRef>>({"gedGsfElectronValueMapsTmp", valueMapPFCandEle_});
0081
0082   valueMapPhoPFCandIso_ = conf_.getParameter<std::string>("valueMapPhoPFblockIso");
0083   valueMapElePFCandIso_ = conf_.getParameter<std::string>("valueMapElePFblockIso");
0084
0085   produces<edm::ValueMap<std::vector<reco::PFCandidateRef>>>(valueMapPhoPFCandIso_);
0086   produces<edm::ValueMap<std::vector<reco::PFCandidateRef>>>(valueMapElePFCandIso_);
0087 }
0088
0089 void ParticleBasedIsoProducer::beginRun(const edm::Run& run, const edm::EventSetup& c) {
0090   thePFBlockBasedIsolation_ = std::make_unique<PFBlockBasedIsolation>();
0091   edm::ParameterSet pfBlockBasedIsolationSetUp = conf_.getParameter<edm::ParameterSet>("pfBlockBasedIsolationSetUp");
0092   thePFBlockBasedIsolation_->setup(pfBlockBasedIsolationSetUp);
0093 }
0094
0095 void ParticleBasedIsoProducer::produce(edm::Event& theEvent, const edm::EventSetup& c) {
0096   auto photonHandle = theEvent.getHandle(photonProducerT_);
0097   auto photonTmpHandle = theEvent.getHandle(photonTmpProducerT_);
0098   auto electronTmpHandle = theEvent.getHandle(electronTmpProducerT_);
0099   auto electronHandle = theEvent.getHandle(electronProducerT_);
0100   auto pfEGCandidateHandle = theEvent.getHandle(pfEgammaCandidates_);
0101   auto pfCandidateHandle = theEvent.getHandle(pfCandidates_);
0102   auto const& pfEGCandToPhotonMap = theEvent.get(valMapPFCandToPhoton_);
0103   auto const& pfEGCandToElectronMap = theEvent.get(valMapPFCandToEle_);
0104
0105   std::vector<std::vector<reco::PFCandidateRef>> pfCandIsoPairVecPho;
0106
0107   ///// Isolation for photons
0108   //  std::cout << " ParticleBasedIsoProducer  photonHandle size " << photonHandle->size() << std::endl;
0109   for (unsigned int lSC = 0; lSC < photonTmpHandle->size(); lSC++) {
0110     reco::PhotonRef phoRef(reco::PhotonRef(photonTmpHandle, lSC));
0111
0112     // loop over the unbiased candidates to retrieve the ref to the unbiased candidate corresponding to this photon
0113     unsigned nObj = pfEGCandidateHandle->size();
0114     reco::PFCandidateRef pfEGCandRef;
0115
0116     std::vector<reco::PFCandidateRef> pfCandIsoPairPho;
0117     for (unsigned int lCand = 0; lCand < nObj; lCand++) {
0118       pfEGCandRef = reco::PFCandidateRef(pfEGCandidateHandle, lCand);
0119       const reco::PhotonRef& myPho = (pfEGCandToPhotonMap)[pfEGCandRef];
0120
0121       if (myPho.isNonnull()) {
0122         //std::cout << "ParticleBasedIsoProducer photons PF SC " << pfEGCandRef->superClusterRef()->energy() << " Photon SC " << myPho->superCluster()->energy() << std::endl;
0123         if (myPho != phoRef)
0124           continue;
0125         //  std::cout << " ParticleBasedIsoProducer photons This is my egammaunbiased guy energy " <<  pfEGCandRef->superClusterRef()->energy() << std::endl;
0126         pfCandIsoPairPho = thePFBlockBasedIsolation_->calculate(myPho->p4(), pfEGCandRef, pfCandidateHandle);
0127
0128         /////// debug
0129         //  for ( std::vector<reco::PFCandidateRef>::const_iterator iPair=pfCandIsoPairPho.begin(); iPair<pfCandIsoPairPho.end(); iPair++) {
0130         // float dR= deltaR(myPho->eta(),  myPho->phi(), (*iPair)->eta(),  (*iPair)->phi() );
0131         // std::cout << " ParticleBasedIsoProducer photons  checking the pfCand bool pair " << (*iPair)->particleId() << " dR " << dR << " pt " <<  (*iPair)->pt() << std::endl;
0132         //  }
0133       }
0134     }
0135
0136     pfCandIsoPairVecPho.push_back(pfCandIsoPairPho);
0137   }
0138
0139   ////////////isolation for electrons
0140   std::vector<std::vector<reco::PFCandidateRef>> pfCandIsoPairVecEle;
0141   //  std::cout << " ParticleBasedIsoProducer  electronHandle size " << electronHandle->size() << std::endl;
0142   for (unsigned int lSC = 0; lSC < electronTmpHandle->size(); lSC++) {
0143     reco::GsfElectronRef eleRef(reco::GsfElectronRef(electronTmpHandle, lSC));
0144
0145     // loop over the unbiased candidates to retrieve the ref to the unbiased candidate corresponding to this electron
0146     unsigned nObj = pfEGCandidateHandle->size();
0147     reco::PFCandidateRef pfEGCandRef;
0148
0149     std::vector<reco::PFCandidateRef> pfCandIsoPairEle;
0150     for (unsigned int lCand = 0; lCand < nObj; lCand++) {
0151       pfEGCandRef = reco::PFCandidateRef(pfEGCandidateHandle, lCand);
0152       const reco::GsfElectronRef& myEle = (pfEGCandToElectronMap)[pfEGCandRef];
0153
0154       if (myEle.isNonnull()) {
0155         //  std::cout << "ParticleBasedIsoProducer Electorns PF SC " << pfEGCandRef->superClusterRef()->energy() << " Electron SC " << myEle->superCluster()->energy() << std::endl;
0156         if (myEle != eleRef)
0157           continue;
0158
0159         //math::XYZVector candidateMomentum(myEle->p4().px(),myEle->p4().py(),myEle->p4().pz());
0160         //math::XYZVector myDir=candidateMomentum.Unit();
0161         //  std::cout << " ParticleBasedIsoProducer  Electrons This is my egammaunbiased guy energy " <<  pfEGCandRef->superClusterRef()->energy()  << std::endl;
0162         //  std::cout << " Ele  direction " << myDir << " eta " << myEle->eta() << " phi " << myEle->phi() << std::endl;
0163         pfCandIsoPairEle = thePFBlockBasedIsolation_->calculate(myEle->p4(), pfEGCandRef, pfCandidateHandle);
0164         /////// debug
0165         //for ( std::vector<reco::PFCandidateRef>::const_iterator iPair=pfCandIsoPairEle.begin(); iPair<pfCandIsoPairEle.end(); iPair++) {
0166         // float dR= deltaR(myEle->eta(),  myEle->phi(), (*iPair)->eta(),  (*iPair)->phi() );
0167         // std::cout << " ParticleBasedIsoProducer Electron  checking the pfCand bool pair " << (*iPair)->particleId() << " dR " << dR << " pt " <<  (*iPair)->pt() << " eta " << (*iPair)->eta() << " phi " << (*iPair)->phi() <<  std::endl;
0168         //  }
0169       }
0170     }
0171
0172     pfCandIsoPairVecEle.push_back(pfCandIsoPairEle);
0173   }
0174
0175   auto phoToPFCandIsoMap_p = std::make_unique<edm::ValueMap<std::vector<reco::PFCandidateRef>>>();
0176   edm::ValueMap<std::vector<reco::PFCandidateRef>>::Filler fillerPhotons(*phoToPFCandIsoMap_p);
0177
0178   //// fill the isolation value map for photons
0179   fillerPhotons.insert(photonHandle, pfCandIsoPairVecPho.begin(), pfCandIsoPairVecPho.end());
0180   fillerPhotons.fill();
0181   theEvent.put(std::move(phoToPFCandIsoMap_p), valueMapPhoPFCandIso_);
0182
0183   auto eleToPFCandIsoMap_p = std::make_unique<edm::ValueMap<std::vector<reco::PFCandidateRef>>>();
0184   edm::ValueMap<std::vector<reco::PFCandidateRef>>::Filler fillerElectrons(*eleToPFCandIsoMap_p);
0185
0186   //// fill the isolation value map for electrons
0187   fillerElectrons.insert(electronHandle, pfCandIsoPairVecEle.begin(), pfCandIsoPairVecEle.end());
0188   fillerElectrons.fill();
0189   theEvent.put(std::move(eleToPFCandIsoMap_p), valueMapElePFCandIso_);
0190 }
0191
0192 void ParticleBasedIsoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0193   // particleBasedIsolation
0194   edm::ParameterSetDescription desc;
0195   desc.add<std::string>("valueMapEleToEG", "");
0196   desc.add<std::string>("valueMapPhoToEG", "valMapPFEgammaCandToPhoton");
0197   desc.add<edm::InputTag>("electronTmpProducer", {"gedGsfElectronsTmp"});
0198   desc.add<edm::InputTag>("pfCandidates", {"particleFlow"});
0199   desc.add<std::string>("valueMapElePFblockIso", "gedGsfElectrons");
0200   desc.add<edm::InputTag>("electronProducer", {"gedGsfElectrons"});
0201   desc.add<edm::InputTag>("photonTmpProducer", {"gedPhotonsTmp"});
0202   desc.add<edm::InputTag>("pfEgammaCandidates", {"particleFlowEGamma"});
0203   {
0204     edm::ParameterSetDescription psd0;
0205     psd0.add<std::string>("ComponentName", "pfBlockBasedIsolation");
0206     psd0.add<double>("coneSize", 9999999999.);
0207     desc.add<edm::ParameterSetDescription>("pfBlockBasedIsolationSetUp", psd0);
0208   }
0209   desc.add<edm::InputTag>("photonProducer", {"gedPhotons"});
0210   desc.add<std::string>("valueMapPhoPFblockIso", "gedPhotons");
0211   descriptions.add("particleBasedIsolation", desc);
0212   // or use the following to generate the label from the module's C++ type
0213   //descriptions.addWithDefaultLabel(desc);
0214 }