NanoAOD/plugins/GenPartIsoProducer.cc

0001 // system include files
0002 #include <memory>
0003
0004 // user include files
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/stream/EDProducer.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/StreamID.h"
0013
0014 #include "DataFormats/JetReco/interface/GenJet.h"
0015 #include "DataFormats/Common/interface/ValueMap.h"
0016
0017 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0018 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0019 #include "DataFormats/PatCandidates/interface/PackedGenParticle.h"
0020 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0021
0022 #include "DataFormats/Math/interface/deltaR.h"
0023
0024 #include <tuple>
0025 #include <string>
0026 #include <vector>
0027 #include <TLorentzVector.h>
0028
0029 using namespace std;
0030
0031 class GenPartIsoProducer : public edm::stream::EDProducer<> {
0032 public:
0033   explicit GenPartIsoProducer(const edm::ParameterSet& iConfig)
0034       : finalGenParticleToken(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genPart"))),
0035         packedGenParticlesToken(
0036             consumes<pat::PackedGenParticleCollection>(iConfig.getParameter<edm::InputTag>("packedGenPart"))),
0037         additionalPdgId_(iConfig.getParameter<int>("additionalPdgId")) {
0038     produces<edm::ValueMap<float>>();
0039   }
0040   ~GenPartIsoProducer() override;
0041
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043
0044 private:
0045   void produce(edm::Event&, const edm::EventSetup&) override;
0046   float computeIso(TLorentzVector thisPart,
0047                    edm::Handle<pat::PackedGenParticleCollection> packedGenParticles,
0048                    std::set<int> gen_fsrset,
0049                    bool skip_leptons);
0050   std::vector<float> Lepts_RelIso;
0051   edm::EDGetTokenT<reco::GenParticleCollection> finalGenParticleToken;
0052   edm::EDGetTokenT<pat::PackedGenParticleCollection> packedGenParticlesToken;
0053   int additionalPdgId_;
0054 };
0055
0056 GenPartIsoProducer::~GenPartIsoProducer() {}
0057
0058 void GenPartIsoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0059   using namespace edm;
0060
0061   auto finalParticles = iEvent.getHandle(finalGenParticleToken);
0062   auto packedGenParticles = iEvent.getHandle(packedGenParticlesToken);
0063
0064   reco::GenParticleCollection::const_iterator genPart;
0065
0066   Lepts_RelIso.clear();
0067
0068   for (genPart = finalParticles->begin(); genPart != finalParticles->end(); genPart++) {
0069     if (abs(genPart->pdgId()) == 11 || abs(genPart->pdgId()) == 13 || abs(genPart->pdgId()) == 15) {
0070       TLorentzVector lep_dressed;
0071       lep_dressed.SetPtEtaPhiE(genPart->pt(), genPart->eta(), genPart->phi(), genPart->energy());
0072       std::set<int> gen_fsrset;
0073       for (size_t k = 0; k < packedGenParticles->size(); k++) {
0074         if ((*packedGenParticles)[k].status() != 1)
0075           continue;
0076         if ((*packedGenParticles)[k].pdgId() != 22)
0077           continue;
0078         double this_dR_lgamma = reco::deltaR(
0079             genPart->eta(), genPart->phi(), (*packedGenParticles)[k].eta(), (*packedGenParticles)[k].phi());
0080         bool idmatch = false;
0081         if ((*packedGenParticles)[k].mother(0)->pdgId() == genPart->pdgId())
0082           idmatch = true;
0083         const reco::Candidate* mother = (*packedGenParticles)[k].mother(0);
0084         for (size_t m = 0; m < mother->numberOfMothers(); m++) {
0085           if ((*packedGenParticles)[k].mother(m)->pdgId() == genPart->pdgId())
0086             idmatch = true;
0087         }
0088         if (!idmatch)
0089           continue;
0090         if (this_dR_lgamma < 0.3) {
0091           gen_fsrset.insert(k);
0092           TLorentzVector gamma;
0093           gamma.SetPtEtaPhiE((*packedGenParticles)[k].pt(),
0094                              (*packedGenParticles)[k].eta(),
0095                              (*packedGenParticles)[k].phi(),
0096                              (*packedGenParticles)[k].energy());
0097           lep_dressed = lep_dressed + gamma;
0098         }
0099       }
0100       float this_GENiso = 0.0;
0101       TLorentzVector thisLep;
0102       thisLep.SetPtEtaPhiM(lep_dressed.Pt(), lep_dressed.Eta(), lep_dressed.Phi(), lep_dressed.M());
0103       this_GENiso = computeIso(thisLep, packedGenParticles, gen_fsrset, true);
0104       Lepts_RelIso.push_back(this_GENiso);
0105     } else if (abs(genPart->pdgId()) == additionalPdgId_) {
0106       float this_GENiso = 0.0;
0107       std::set<int> gen_fsrset_nolep;
0108       TLorentzVector thisPart;
0109       thisPart.SetPtEtaPhiE(genPart->pt(), genPart->eta(), genPart->phi(), genPart->energy());
0110       this_GENiso = computeIso(thisPart, packedGenParticles, gen_fsrset_nolep, false);
0111       Lepts_RelIso.push_back(this_GENiso);
0112     } else {
0113       float this_GENiso = 0.0;
0114       Lepts_RelIso.push_back(this_GENiso);
0115     }
0116   }
0117
0118   auto isoV = std::make_unique<edm::ValueMap<float>>();
0119   edm::ValueMap<float>::Filler fillerIsoMap(*isoV);
0120   fillerIsoMap.insert(finalParticles, Lepts_RelIso.begin(), Lepts_RelIso.end());
0121   fillerIsoMap.fill();
0122   iEvent.put(std::move(isoV));
0123 }
0124
0125 float GenPartIsoProducer::computeIso(TLorentzVector thisPart,
0126                                      edm::Handle<pat::PackedGenParticleCollection> packedGenParticles,
0127                                      std::set<int> gen_fsrset,
0128                                      bool skip_leptons) {
0129   double this_GENiso = 0.0;
0130   for (size_t k = 0; k < packedGenParticles->size(); k++) {
0131     if ((*packedGenParticles)[k].status() != 1)
0132       continue;
0133     if (abs((*packedGenParticles)[k].pdgId()) == 12 || abs((*packedGenParticles)[k].pdgId()) == 14 ||
0134         abs((*packedGenParticles)[k].pdgId()) == 16)
0135       continue;
0136     if (reco::deltaR(thisPart.Eta(), thisPart.Phi(), (*packedGenParticles)[k].eta(), (*packedGenParticles)[k].phi()) <
0137         0.001)
0138       continue;
0139     if (skip_leptons == true) {
0140       if ((abs((*packedGenParticles)[k].pdgId()) == 11 || abs((*packedGenParticles)[k].pdgId()) == 13))
0141         continue;
0142       if (gen_fsrset.find(k) != gen_fsrset.end())
0143         continue;
0144     }
0145     double this_dRvL_nolep =
0146         reco::deltaR(thisPart.Eta(), thisPart.Phi(), (*packedGenParticles)[k].eta(), (*packedGenParticles)[k].phi());
0147     if (this_dRvL_nolep < 0.3) {
0148       this_GENiso = this_GENiso + (*packedGenParticles)[k].pt();
0149     }
0150   }
0151   this_GENiso = this_GENiso / thisPart.Pt();
0152   return this_GENiso;
0153 }
0154
0155 void GenPartIsoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0156   // Description of external inputs requered by this module
0157   // genPart: collection of gen particles for which to compute Iso
0158   // packedGenPart: collection of particles to be used for leptons dressing
0159   // additionalPdgId: additional particle (besides leptons) for which Iso is computed
0160   edm::ParameterSetDescription desc;
0161   desc.add<edm::InputTag>("genPart")->setComment("input physics object collection");
0162   desc.add<edm::InputTag>("packedGenPart")->setComment("input stable hadrons collection");
0163   desc.add<int>("additionalPdgId")->setComment("additional pdgId for Iso computation");
0164   descriptions.addDefault(desc);
0165 }
0166
0167 //define this as a plug-in
0168 DEFINE_FWK_MODULE(GenPartIsoProducer);