diff --git a/PhysicsTools/NanoAOD/plugins/GenPartIsoProducer.cc b/PhysicsTools/NanoAOD/plugins/GenPartIsoProducer.cc new file mode 100644 index 0000000000000..5ba16daa947ee --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/GenPartIsoProducer.cc @@ -0,0 +1,168 @@ +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/JetReco/interface/GenJet.h" +#include "DataFormats/Common/interface/ValueMap.h" + +#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h" +#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" +#include "DataFormats/PatCandidates/interface/PackedGenParticle.h" +#include "DataFormats/PatCandidates/interface/CompositeCandidate.h" + +#include "DataFormats/Math/interface/deltaR.h" + +#include +#include +#include +#include + +using namespace std; + +class GenPartIsoProducer : public edm::stream::EDProducer<> { +public: + explicit GenPartIsoProducer(const edm::ParameterSet& iConfig) + : finalGenParticleToken(consumes(iConfig.getParameter("genPart"))), + packedGenParticlesToken( + consumes(iConfig.getParameter("packedGenPart"))), + additionalPdgId_(iConfig.getParameter("additionalPdgId")) { + produces>(); + } + ~GenPartIsoProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + float computeIso(TLorentzVector thisPart, + edm::Handle packedGenParticles, + std::set gen_fsrset, + bool skip_leptons); + std::vector Lepts_RelIso; + edm::EDGetTokenT finalGenParticleToken; + edm::EDGetTokenT packedGenParticlesToken; + int additionalPdgId_; +}; + +GenPartIsoProducer::~GenPartIsoProducer() {} + +void GenPartIsoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + using namespace edm; + + auto finalParticles = iEvent.getHandle(finalGenParticleToken); + auto packedGenParticles = iEvent.getHandle(packedGenParticlesToken); + + reco::GenParticleCollection::const_iterator genPart; + + Lepts_RelIso.clear(); + + for (genPart = finalParticles->begin(); genPart != finalParticles->end(); genPart++) { + if (abs(genPart->pdgId()) == 11 || abs(genPart->pdgId()) == 13 || abs(genPart->pdgId()) == 15) { + TLorentzVector lep_dressed; + lep_dressed.SetPtEtaPhiE(genPart->pt(), genPart->eta(), genPart->phi(), genPart->energy()); + std::set gen_fsrset; + for (size_t k = 0; k < packedGenParticles->size(); k++) { + if ((*packedGenParticles)[k].status() != 1) + continue; + if ((*packedGenParticles)[k].pdgId() != 22) + continue; + double this_dR_lgamma = reco::deltaR( + genPart->eta(), genPart->phi(), (*packedGenParticles)[k].eta(), (*packedGenParticles)[k].phi()); + bool idmatch = false; + if ((*packedGenParticles)[k].mother(0)->pdgId() == genPart->pdgId()) + idmatch = true; + const reco::Candidate* mother = (*packedGenParticles)[k].mother(0); + for (size_t m = 0; m < mother->numberOfMothers(); m++) { + if ((*packedGenParticles)[k].mother(m)->pdgId() == genPart->pdgId()) + idmatch = true; + } + if (!idmatch) + continue; + if (this_dR_lgamma < 0.3) { + gen_fsrset.insert(k); + TLorentzVector gamma; + gamma.SetPtEtaPhiE((*packedGenParticles)[k].pt(), + (*packedGenParticles)[k].eta(), + (*packedGenParticles)[k].phi(), + (*packedGenParticles)[k].energy()); + lep_dressed = lep_dressed + gamma; + } + } + float this_GENiso = 0.0; + TLorentzVector thisLep; + thisLep.SetPtEtaPhiM(lep_dressed.Pt(), lep_dressed.Eta(), lep_dressed.Phi(), lep_dressed.M()); + this_GENiso = computeIso(thisLep, packedGenParticles, gen_fsrset, true); + Lepts_RelIso.push_back(this_GENiso); + } else if (abs(genPart->pdgId()) == additionalPdgId_) { + float this_GENiso = 0.0; + std::set gen_fsrset_nolep; + TLorentzVector thisPart; + thisPart.SetPtEtaPhiE(genPart->pt(), genPart->eta(), genPart->phi(), genPart->energy()); + this_GENiso = computeIso(thisPart, packedGenParticles, gen_fsrset_nolep, false); + Lepts_RelIso.push_back(this_GENiso); + } else { + float this_GENiso = 0.0; + Lepts_RelIso.push_back(this_GENiso); + } + } + + auto isoV = std::make_unique>(); + edm::ValueMap::Filler fillerIsoMap(*isoV); + fillerIsoMap.insert(finalParticles, Lepts_RelIso.begin(), Lepts_RelIso.end()); + fillerIsoMap.fill(); + iEvent.put(std::move(isoV)); +} + +float GenPartIsoProducer::computeIso(TLorentzVector thisPart, + edm::Handle packedGenParticles, + std::set gen_fsrset, + bool skip_leptons) { + double this_GENiso = 0.0; + for (size_t k = 0; k < packedGenParticles->size(); k++) { + if ((*packedGenParticles)[k].status() != 1) + continue; + if (abs((*packedGenParticles)[k].pdgId()) == 12 || abs((*packedGenParticles)[k].pdgId()) == 14 || + abs((*packedGenParticles)[k].pdgId()) == 16) + continue; + if (reco::deltaR(thisPart.Eta(), thisPart.Phi(), (*packedGenParticles)[k].eta(), (*packedGenParticles)[k].phi()) < + 0.001) + continue; + if (skip_leptons == true) { + if ((abs((*packedGenParticles)[k].pdgId()) == 11 || abs((*packedGenParticles)[k].pdgId()) == 13)) + continue; + if (gen_fsrset.find(k) != gen_fsrset.end()) + continue; + } + double this_dRvL_nolep = + reco::deltaR(thisPart.Eta(), thisPart.Phi(), (*packedGenParticles)[k].eta(), (*packedGenParticles)[k].phi()); + if (this_dRvL_nolep < 0.3) { + this_GENiso = this_GENiso + (*packedGenParticles)[k].pt(); + } + } + this_GENiso = this_GENiso / thisPart.Pt(); + return this_GENiso; +} + +void GenPartIsoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + // Description of external inputs requered by this module + // genPart: collection of gen particles for which to compute Iso + // packedGenPart: collection of particles to be used for leptons dressing + // additionalPdgId: additional particle (besides leptons) for which Iso is computed + edm::ParameterSetDescription desc; + desc.add("genPart")->setComment("input physics object collection"); + desc.add("packedGenPart")->setComment("input stable hadrons collection"); + desc.add("additionalPdgId")->setComment("additional pdgId for Iso computation"); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(GenPartIsoProducer); diff --git a/PhysicsTools/NanoAOD/python/genparticles_cff.py b/PhysicsTools/NanoAOD/python/genparticles_cff.py index d5191f14cf290..33e09275d7c97 100644 --- a/PhysicsTools/NanoAOD/python/genparticles_cff.py +++ b/PhysicsTools/NanoAOD/python/genparticles_cff.py @@ -8,22 +8,22 @@ finalGenParticles = cms.EDProducer("GenParticlePruner", src = cms.InputTag("prunedGenParticles"), select = cms.vstring( - "drop *", + "drop *", "keep++ abs(pdgId) == 15 & (pt > 15 || isPromptDecayed() )",# keep full tau decay chain for some taus - #"drop status==1 & pt < 1", #drop soft stable particle in tau decay + #"drop status==1 & pt < 1", #drop soft stable particle in tau decay "keep+ abs(pdgId) == 15 ", # keep first gen decay product for all tau "+keep pdgId == 22 && status == 1 && (pt > 10 || isPromptFinalState())", # keep gamma above 10 GeV (or all prompt) and its first parent - "+keep abs(pdgId) == 11 || abs(pdgId) == 13 || abs(pdgId) == 15", #keep leptons, with at most one mother back in the history - "drop abs(pdgId)= 2212 && abs(pz) > 1000", #drop LHC protons accidentally added by previous keeps + "+keep abs(pdgId) == 11 || abs(pdgId) == 13 || abs(pdgId) == 15", #keep leptons, with at most one mother back in the history + "drop abs(pdgId)= 2212 && abs(pz) > 1000", #drop LHC protons accidentally added by previous keeps "keep (400 < abs(pdgId) < 600) || (4000 < abs(pdgId) < 6000)", #keep all B and C hadrons "keep abs(pdgId) == 12 || abs(pdgId) == 14 || abs(pdgId) == 16", # keep neutrinos - "keep status == 3 || (status > 20 && status < 30)", #keep matrix element summary + "keep status == 3 || (status > 20 && status < 30)", #keep matrix element summary "keep isHardProcess() || fromHardProcessDecayed() || fromHardProcessFinalState() || (statusFlags().fromHardProcess() && statusFlags().isLastCopy())", #keep event summary based on status flags - "keep (status > 70 && status < 80 && pt > 15) ", # keep high pt partons right before hadronization + "keep (status > 70 && status < 80 && pt > 15) ", # keep high pt partons right before hadronization "keep abs(pdgId) == 23 || abs(pdgId) == 24 || abs(pdgId) == 25 || abs(pdgId) == 37 ", # keep VIP(articles)s #"keep abs(pdgId) == 310 && abs(eta) < 2.5 && pt > 1 ", # keep K0 "keep (1000001 <= abs(pdgId) <= 1000039 ) || ( 2000001 <= abs(pdgId) <= 2000015)", #keep SUSY fiction particles - ) + ) ) @@ -33,6 +33,9 @@ src = cms.InputTag("finalGenParticles"), name= cms.string("GenPart"), doc = cms.string("interesting gen particles "), + externalVariables = cms.PSet( + iso = ExtVar(cms.InputTag("genIso"), float, precision=8, doc="Isolation for leptons"), + ), variables = cms.PSet( pt = Var("pt", float, precision=8), phi = Var("phi", float,precision=8), @@ -77,6 +80,11 @@ ) ) -genParticleTask = cms.Task(finalGenParticles) -genParticleTablesTask = cms.Task(genParticleTable) +genIso = cms.EDProducer("GenPartIsoProducer", + genPart=cms.InputTag("finalGenParticles"), + packedGenPart=cms.InputTag("packedGenParticles"), + additionalPdgId=cms.int32(22), + ) +genParticleTask = cms.Task(finalGenParticles, genIso) +genParticleTablesTask = cms.Task(genParticleTable) diff --git a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py index 6634327de7b58..066333e589ca9 100644 --- a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py +++ b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py @@ -351,6 +351,7 @@ Plot1D('pdgId', 'pdgId', 20, -6000, 6000, 'PDG id'), Plot1D('phi', 'phi', 20, -3.14159, 3.14159, 'phi'), Plot1D('pt', 'pt', 20, 0, 200, 'pt'), + Plot1D('iso', 'iso', 20, 0, 200, 'iso'), Plot1D('status', 'status', 20, 0, 100, 'Particle status. 1=stable'), Plot1D('statusFlags', 'statusFlags', 15, 0, 15, 'gen status flags stored bitwise, bits are: 0 : isPrompt, 1 : isDecayedLeptonHadron, 2 : isTauDecayProduct, 3 : isPromptTauDecayProduct, 4 : isDirectTauDecayProduct, 5 : isDirectPromptTauDecayProduct, 6 : isDirectHadronDecayProduct, 7 : isHardProcess, 8 : fromHardProcess, 9 : isHardProcessTauDecayProduct, 10 : isDirectHardProcessTauDecayProduct, 11 : fromHardProcessBeforeFSR, 12 : isFirstCopy, 13 : isLastCopy, 14 : isLastCopyBeforeFSR, ', bitset=True), ) diff --git a/PhysicsTools/NanoAOD/python/nanogen_cff.py b/PhysicsTools/NanoAOD/python/nanogen_cff.py index 8c4915b78417e..45a8059ee40f2 100644 --- a/PhysicsTools/NanoAOD/python/nanogen_cff.py +++ b/PhysicsTools/NanoAOD/python/nanogen_cff.py @@ -26,6 +26,7 @@ genJetAK8FlavourTable+ cms.Sequence(genTauTask)+ genTable+ + genIso+ genFilterTable+ cms.Sequence(genParticleTablesTask)+ cms.Sequence(genVertexTablesTask)+ @@ -97,6 +98,8 @@ def customizeNanoGEN(process): process.nanogenSequence.remove(process.genParticles2HepMCHiggsVtx) process.nanogenSequence.remove(process.genParticles2HepMC) process.nanogenSequence.remove(process.mergedGenParticles) + process.nanogenSequence.remove(process.genIso) + delattr(process.genParticleTable.externalVariables,"iso") nanoGenCommonCustomize(process) return process