diff --git a/GeneratorInterface/GenFilters/plugins/AJJGenJetFilter.cc b/GeneratorInterface/GenFilters/plugins/AJJGenJetFilter.cc new file mode 100644 index 0000000000000..0287a1d87e3e7 --- /dev/null +++ b/GeneratorInterface/GenFilters/plugins/AJJGenJetFilter.cc @@ -0,0 +1,241 @@ +// -*- C++ -*- +// +// Package: GeneratorInterface/GenFilters +// Class: AJJGenJetFilter +// +/* + + Description: A filter to select events with one photon and 2 jets. + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: Hamed Bakhshian +// Created: Wed Oct 06 2021 +// +// + +// CMSSW include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDFilter.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" + +#include "DataFormats/JetReco/interface/GenJetCollection.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" + +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "DataFormats/Math/interface/deltaPhi.h" +#include "DataFormats/Math/interface/deltaR.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include + +// C++ include files +#include +#include +#include +#include + +using namespace edm; +using namespace std; +// +// class declaration +// + +class AJJGenJetFilter : public edm::global::EDFilter<> { +public: + explicit AJJGenJetFilter(const edm::ParameterSet& pset); + ~AJJGenJetFilter() override; + + static void fillDescriptions(edm::ConfigurationDescriptions&); + bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + +private: + // ----------memeber function---------------------- + const std::vector filterGenJets(const vector* jets) const; + const std::vector filterGenLeptons(const std::vector* particles) const; + const std::vector filterGenPhotons(const std::vector* particles) const; + + //************************** + // Private Member data ***** +private: + // Dijet cut + const double ptMin; + const double etaMin; + const double etaMax; + const double deltaRJetLep; + + const double minDeltaEta; + const double maxDeltaEta; + const double mininvmass; + + const double maxPhotonEta; + const double minPhotonPt; + const double maxPhotonPt; + + // Input tags + edm::EDGetTokenT m_GenJetCollection; + edm::EDGetTokenT m_GenParticleCollection; +}; + +AJJGenJetFilter::AJJGenJetFilter(const edm::ParameterSet& iConfig) + : ptMin(iConfig.getParameter("minPt")), + etaMin(iConfig.getParameter("minEta")), + etaMax(iConfig.getParameter("maxEta")), + deltaRJetLep(iConfig.getParameter("deltaRJetLep")), + + minDeltaEta(iConfig.getParameter("minDeltaEta")), + maxDeltaEta(iConfig.getParameter("maxDeltaEta")), + mininvmass(iConfig.getParameter("MinInvMass")), + + maxPhotonEta(iConfig.getParameter("maxPhotonEta")), + minPhotonPt(iConfig.getParameter("minPhotonPt")), + maxPhotonPt(iConfig.getParameter("maxPhotonPt")) { + m_GenParticleCollection = consumes(iConfig.getParameter("genParticles")); + + if (ptMin > 0) { + m_GenJetCollection = consumes(iConfig.getParameter("GenJetCollection")); + } + edm::LogInfo("AJJGenJetFilter") << "Parameters:" + << "jPtMin:" << ptMin << ",jEta:" << etaMin << "--" << etaMax << ",minDR(j,lep)" + << deltaRJetLep << ",deltaEta(j1,j2)" << minDeltaEta << "--" << maxDeltaEta + << "m(j1,j2) < " << mininvmass << "PhotonPt" << minPhotonPt << "--" << maxPhotonPt + << "PhotonEta" << maxPhotonEta; +} + +AJJGenJetFilter::~AJJGenJetFilter() {} + +const vector AJJGenJetFilter::filterGenLeptons( + const vector* particles) const { + vector out; + + for (const auto& p : *particles) { + int absPdgId = std::abs(p.pdgId()); + + if (((absPdgId == 11) || (absPdgId == 13) || (absPdgId == 15)) && p.isHardProcess()) { + out.push_back(&p); + } + } + return out; +} + +const vector AJJGenJetFilter::filterGenPhotons( + const vector* particles) const { + vector out; + + for (const auto& p : *particles) { + int absPdgId = std::abs(p.pdgId()); + + if ((absPdgId == 22) && p.isHardProcess()) { + if (std::abs(p.eta()) < maxPhotonEta && p.pt() > minPhotonPt && p.pt() <= maxPhotonPt) { + out.push_back(&p); + } else { + edm::LogInfo("AJJPhoton") << "photon rejected, pt:" << p.pt() << " , eta:" << p.eta(); + } + } + } + + return out; +} + +const vector AJJGenJetFilter::filterGenJets(const vector* jets) const { + vector out; + + for (auto const& j : *jets) { + if (j.p4().pt() > ptMin && j.p4().eta() > etaMin && j.p4().eta() < etaMax) { + out.push_back(&j); + } else { + edm::LogInfo("AJJJets") << "Jet rejected, pt:" << j.p4().pt() << " eta:" << j.p4().eta(); + } + } + + return out; +} + +// ------------ method called to skim the data ------------ +bool AJJGenJetFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const { + using namespace edm; + // Getting filtered generator jets + + //Handle genParticelesCollection; + //iEvent.getByToken(m_GenParticleCollection, genParticelesCollection); + //const vector* genParticles = genParticelesCollection.product(); + auto const& genParticles = iEvent.get(m_GenParticleCollection); + vector filGenLep = filterGenLeptons(&genParticles); + + vector filGenPhotons = filterGenPhotons(&genParticles); + + if (filGenPhotons.size() != 1) { + edm::LogInfo("AJJPhoton") << "Events rejected, number of photons:" << filGenPhotons.size(); + return false; + } + + if (ptMin < 0) + return true; + + //Handle > handleGenJets; + //iEvent.getByToken(m_GenJetCollection, handleGenJets); + //const vector* genJets = handleGenJets.product(); + auto const& genJets = iEvent.get(m_GenJetCollection); + + vector filGenJets = filterGenJets(&genJets); + // Getting p4 of jet with no lepton + vector genJetsWithoutLeptonsP4; + unsigned int nGoodJets = 0; + for (auto const& j : filGenJets) { + bool cleanJet = true; + const math::XYZTLorentzVector& p4J = j->p4(); + for (auto const& p : filGenLep) + if (reco::deltaR2(p->p4(), p4J) < deltaRJetLep * deltaRJetLep) + cleanJet = false; + + if (cleanJet) { + if (genJetsWithoutLeptonsP4.size() < 2) + genJetsWithoutLeptonsP4.push_back(p4J); + nGoodJets++; + } + } + + //If we do not find at least 2 jets veto the event + if (nGoodJets < 2) { + edm::LogInfo("AJJJets") << "Events rejected, number of jets:" << nGoodJets; + return false; + } + + double dEta = std::abs(genJetsWithoutLeptonsP4[0].eta() - genJetsWithoutLeptonsP4[1].eta()); + float invMassLeadingJet = (genJetsWithoutLeptonsP4[0] + genJetsWithoutLeptonsP4[1]).M(); + + if (dEta >= minDeltaEta && dEta <= maxDeltaEta && invMassLeadingJet > mininvmass) { + return true; + } + + edm::LogInfo("AJJJets") << "Events rejected, dEta:" << dEta << ", mjj:" << invMassLeadingJet; + return false; +} + +void AJJGenJetFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("GenJetCollection", edm::InputTag("ak4GenJetsNoNu")); + desc.add("genParticles", edm::InputTag("genParticles")); + desc.add("minPt", -1.0)->setComment("If this is negative, no cut on jets is applied"); + desc.add("minEta", -5.0); + desc.add("maxEta", 5.0); + desc.add("deltaRJetLep", 0.); + desc.add("minDeltaEta", 0.0); + desc.add("maxDeltaEta", 9999.0); + desc.add("MinInvMass", 0.0); + desc.add("maxPhotonEta", 5); + desc.add("minPhotonPt", 50); + desc.add("maxPhotonPt", 10000); + descriptions.addDefault(desc); +} +//define this as a plug-in +DEFINE_FWK_MODULE(AJJGenJetFilter); diff --git a/GeneratorInterface/GenFilters/python/AJJGenJetFilter_cff.py b/GeneratorInterface/GenFilters/python/AJJGenJetFilter_cff.py new file mode 100644 index 0000000000000..ece4603364f02 --- /dev/null +++ b/GeneratorInterface/GenFilters/python/AJJGenJetFilter_cff.py @@ -0,0 +1,12 @@ +import FWCore.ParameterSet.Config as cms + +from GeneratorInterface.GenFilters.AJJGenJetFilter_cfi import * +from RecoJets.Configuration.GenJetParticles_cff import * +from RecoJets.Configuration.RecoGenJets_cff import * + +vjjGenJetFilterSeq = cms.Sequence( + genParticlesForJetsNoNu* + ak4GenJetsNoNu* + vjjGenJetFilterPhotonInBarrelMjj300 +) + diff --git a/GeneratorInterface/GenFilters/python/AJJGenJetFilter_cfi.py b/GeneratorInterface/GenFilters/python/AJJGenJetFilter_cfi.py new file mode 100644 index 0000000000000..e24b523fc94ac --- /dev/null +++ b/GeneratorInterface/GenFilters/python/AJJGenJetFilter_cfi.py @@ -0,0 +1,44 @@ +import FWCore.ParameterSet.Config as cms + +ajjGenJetFilterPhotonInBarrelMjj300 = cms.EDFilter("AJJGenJetFilter", + GenJetCollection = cms.InputTag('ak4GenJetsNoNu'), + genParticles = cms.InputTag("genParticles"), + + #following cuts are applied to select jets + minPt = cms.untracked.double( 40 ), + minEta = cms.untracked.double( -4.5 ), + maxEta = cms.untracked.double( 4.5 ), + deltaRJetLep = cms.untracked.double( 0.3 ), + + #following cuts are applied on the first two leading jets + minDeltaEta = cms.untracked.double( 3.0 ), + maxDeltaEta = cms.untracked.double( 999.0 ), + MinInvMass = cms.untracked.double( 300 ), + + #the cut on the photon eta + maxPhotonEta = cms.untracked.double( 1.48 ), + minPhotonPt = cms.untracked.double( 50 ), + maxPhotonPt = cms.untracked.double( 10000 ) +) + +ajjGenJetFilterPhoton = cms.EDFilter("AJJGenJetFilter", + GenJetCollection = cms.InputTag('ak4GenJetsNoNu'), + genParticles = cms.InputTag("genParticles"), + + #following cuts are applied to select jets + #if minPt is negative, no criteri on jets (including njets and delta_eta and invmasss) is applied + minPt = cms.untracked.double( -1 ), + minEta = cms.untracked.double( -4.5 ), + maxEta = cms.untracked.double( 4.5 ), + deltaRJetLep = cms.untracked.double( 0.3 ), + + #following cuts are applied on the first two leading jets + minDeltaEta = cms.untracked.double( 3.0 ), + maxDeltaEta = cms.untracked.double( 999.0 ), + MinInvMass = cms.untracked.double( 300 ), + + #the cut on the photon eta + maxPhotonEta = cms.untracked.double( 1.48 ), + minPhotonPt = cms.untracked.double( 50 ), + maxPhotonPt = cms.untracked.double( 10000 ) +)