Skip to content

Commit

Permalink
Merge pull request #44654 from namapane/dev_geniso_140X
Browse files Browse the repository at this point in the history
[14_0_X] New plugin for gen-part isolation.
  • Loading branch information
cmsbuild authored Apr 12, 2024
2 parents d39f442 + 4ee621c commit 050099f
Show file tree
Hide file tree
Showing 4 changed files with 189 additions and 9 deletions.
168 changes: 168 additions & 0 deletions PhysicsTools/NanoAOD/plugins/GenPartIsoProducer.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,168 @@
// system include files
#include <memory>

// 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 <tuple>
#include <string>
#include <vector>
#include <TLorentzVector.h>

using namespace std;

class GenPartIsoProducer : public edm::stream::EDProducer<> {
public:
explicit GenPartIsoProducer(const edm::ParameterSet& iConfig)
: finalGenParticleToken(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genPart"))),
packedGenParticlesToken(
consumes<pat::PackedGenParticleCollection>(iConfig.getParameter<edm::InputTag>("packedGenPart"))),
additionalPdgId_(iConfig.getParameter<int>("additionalPdgId")) {
produces<edm::ValueMap<float>>();
}
~GenPartIsoProducer() override;

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
void produce(edm::Event&, const edm::EventSetup&) override;
float computeIso(TLorentzVector thisPart,
edm::Handle<pat::PackedGenParticleCollection> packedGenParticles,
std::set<int> gen_fsrset,
bool skip_leptons);
std::vector<float> Lepts_RelIso;
edm::EDGetTokenT<reco::GenParticleCollection> finalGenParticleToken;
edm::EDGetTokenT<pat::PackedGenParticleCollection> 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<int> 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<int> 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<float>>();
edm::ValueMap<float>::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<pat::PackedGenParticleCollection> packedGenParticles,
std::set<int> 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<edm::InputTag>("genPart")->setComment("input physics object collection");
desc.add<edm::InputTag>("packedGenPart")->setComment("input stable hadrons collection");
desc.add<int>("additionalPdgId")->setComment("additional pdgId for Iso computation");
descriptions.addDefault(desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(GenPartIsoProducer);
26 changes: 17 additions & 9 deletions PhysicsTools/NanoAOD/python/genparticles_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
)
)


Expand All @@ -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),
Expand Down Expand Up @@ -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)
1 change: 1 addition & 0 deletions PhysicsTools/NanoAOD/python/nanoDQM_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -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),
)
Expand Down
3 changes: 3 additions & 0 deletions PhysicsTools/NanoAOD/python/nanogen_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
genJetAK8FlavourTable+
cms.Sequence(genTauTask)+
genTable+
genIso+
genFilterTable+
cms.Sequence(genParticleTablesTask)+
cms.Sequence(genVertexTablesTask)+
Expand Down Expand Up @@ -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

Expand Down

0 comments on commit 050099f

Please sign in to comment.