GenFilters/plugins/AJJGenJetFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    GeneratorInterface/GenFilters
0004 // Class:      AJJGenJetFilter
0005 //
0006 /*
0007
0008  Description: A filter to select events with one photon and 2 jets.
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Hamed Bakhshian
0015 //         Created:  Wed Oct 06 2021
0016 //
0017 //
0018
0019 // CMSSW include files
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/global/EDFilter.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026
0027 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0028 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0029
0030 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0031 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0032
0033 #include "FWCore/Utilities/interface/EDGetToken.h"
0034 #include "FWCore/Utilities/interface/InputTag.h"
0035 #include "DataFormats/Math/interface/deltaPhi.h"
0036 #include "DataFormats/Math/interface/deltaR.h"
0037
0038 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0039 #include <HepMC/GenVertex.h>
0040
0041 // C++ include files
0042 #include <memory>
0043 #include <map>
0044 #include <vector>
0045 #include <iostream>
0046
0047 using namespace edm;
0048 using namespace std;
0049 //
0050 // class declaration
0051 //
0052
0053 class AJJGenJetFilter : public edm::global::EDFilter<> {
0054 public:
0055   explicit AJJGenJetFilter(const edm::ParameterSet& pset);
0056   ~AJJGenJetFilter() override;
0057
0058   static void fillDescriptions(edm::ConfigurationDescriptions&);
0059   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0060
0061 private:
0062   // ----------memeber function----------------------
0063   const std::vector<const reco::GenJet*> filterGenJets(const vector<reco::GenJet>* jets) const;
0064   const std::vector<const reco::GenParticle*> filterGenLeptons(const std::vector<reco::GenParticle>* particles) const;
0065   const std::vector<const reco::GenParticle*> filterGenPhotons(const std::vector<reco::GenParticle>* particles) const;
0066
0067   //**************************
0068   // Private Member data *****
0069 private:
0070   // Dijet cut
0071   const double ptMin;
0072   const double etaMin;
0073   const double etaMax;
0074   const double deltaRJetLep;
0075
0076   const double minDeltaEta;
0077   const double maxDeltaEta;
0078   const double mininvmass;
0079
0080   const double maxPhotonEta;
0081   const double minPhotonPt;
0082   const double maxPhotonPt;
0083
0084   // Input tags
0085   edm::EDGetTokenT<reco::GenJetCollection> m_GenJetCollection;
0086   edm::EDGetTokenT<reco::GenParticleCollection> m_GenParticleCollection;
0087 };
0088
0089 AJJGenJetFilter::AJJGenJetFilter(const edm::ParameterSet& iConfig)
0090     : ptMin(iConfig.getParameter<double>("minPt")),
0091       etaMin(iConfig.getParameter<double>("minEta")),
0092       etaMax(iConfig.getParameter<double>("maxEta")),
0093       deltaRJetLep(iConfig.getParameter<double>("deltaRJetLep")),
0094
0095       minDeltaEta(iConfig.getParameter<double>("minDeltaEta")),
0096       maxDeltaEta(iConfig.getParameter<double>("maxDeltaEta")),
0097       mininvmass(iConfig.getParameter<double>("MinInvMass")),
0098
0099       maxPhotonEta(iConfig.getParameter<double>("maxPhotonEta")),
0100       minPhotonPt(iConfig.getParameter<double>("minPhotonPt")),
0101       maxPhotonPt(iConfig.getParameter<double>("maxPhotonPt")) {
0102   m_GenParticleCollection = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticles"));
0103
0104   if (ptMin > 0) {
0105     m_GenJetCollection = consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("GenJetCollection"));
0106   }
0107   edm::LogInfo("AJJGenJetFilter") << "Parameters:"
0108                                   << "jPtMin:" << ptMin << ",jEta:" << etaMin << "--" << etaMax << ",minDR(j,lep)"
0109                                   << deltaRJetLep << ",deltaEta(j1,j2)" << minDeltaEta << "--" << maxDeltaEta
0110                                   << "m(j1,j2) < " << mininvmass << "PhotonPt" << minPhotonPt << "--" << maxPhotonPt
0111                                   << "PhotonEta" << maxPhotonEta;
0112 }
0113
0114 AJJGenJetFilter::~AJJGenJetFilter() {}
0115
0116 const vector<const reco::GenParticle*> AJJGenJetFilter::filterGenLeptons(
0117     const vector<reco::GenParticle>* particles) const {
0118   vector<const reco::GenParticle*> out;
0119
0120   for (const auto& p : *particles) {
0121     int absPdgId = std::abs(p.pdgId());
0122
0123     if (((absPdgId == 11) || (absPdgId == 13) || (absPdgId == 15)) && p.isHardProcess()) {
0124       out.push_back(&p);
0125     }
0126   }
0127   return out;
0128 }
0129
0130 const vector<const reco::GenParticle*> AJJGenJetFilter::filterGenPhotons(
0131     const vector<reco::GenParticle>* particles) const {
0132   vector<const reco::GenParticle*> out;
0133
0134   for (const auto& p : *particles) {
0135     int absPdgId = std::abs(p.pdgId());
0136
0137     if ((absPdgId == 22) && p.isHardProcess()) {
0138       if (std::abs(p.eta()) < maxPhotonEta && p.pt() > minPhotonPt && p.pt() <= maxPhotonPt) {
0139         out.push_back(&p);
0140       } else {
0141         edm::LogInfo("AJJPhoton") << "photon rejected, pt:" << p.pt() << " , eta:" << p.eta();
0142       }
0143     }
0144   }
0145
0146   return out;
0147 }
0148
0149 const vector<const reco::GenJet*> AJJGenJetFilter::filterGenJets(const vector<reco::GenJet>* jets) const {
0150   vector<const reco::GenJet*> out;
0151
0152   for (auto const& j : *jets) {
0153     if (j.p4().pt() > ptMin && j.p4().eta() > etaMin && j.p4().eta() < etaMax) {
0154       out.push_back(&j);
0155     } else {
0156       edm::LogInfo("AJJJets") << "Jet rejected, pt:" << j.p4().pt() << " eta:" << j.p4().eta();
0157     }
0158   }
0159
0160   return out;
0161 }
0162
0163 // ------------ method called to skim the data  ------------
0164 bool AJJGenJetFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0165   using namespace edm;
0166   // Getting filtered generator jets
0167
0168   auto const& genParticles = iEvent.get(m_GenParticleCollection);
0169
0170   vector<const reco::GenParticle*> filGenPhotons = filterGenPhotons(&genParticles);
0171
0172   if (filGenPhotons.size() != 1) {
0173     edm::LogInfo("AJJPhoton") << "Events rejected, number of photons:" << filGenPhotons.size();
0174     return false;
0175   }
0176
0177   if (ptMin < 0)
0178     return true;
0179
0180   vector<const reco::GenParticle*> filGenLep = filterGenLeptons(&genParticles);
0181   auto const& genJets = iEvent.get(m_GenJetCollection);
0182
0183   vector<const reco::GenJet*> filGenJets = filterGenJets(&genJets);
0184   // Getting p4 of jet with no lepton
0185   vector<math::XYZTLorentzVector> genJetsWithoutLeptonsP4;
0186   unsigned int nGoodJets = 0;
0187   for (auto const& j : filGenJets) {
0188     bool cleanJet = true;
0189     const math::XYZTLorentzVector& p4J = j->p4();
0190     for (auto const& p : filGenLep)
0191       if (reco::deltaR2(p->p4(), p4J) < deltaRJetLep * deltaRJetLep) {
0192         cleanJet = false;
0193         break;
0194       }
0195     if (cleanJet) {
0196       if (genJetsWithoutLeptonsP4.size() < 2)
0197         genJetsWithoutLeptonsP4.push_back(p4J);
0198       nGoodJets++;
0199     }
0200   }
0201
0202   //If we do not find at least 2 jets veto the event
0203   if (nGoodJets < 2) {
0204     LogDebug("AJJJets") << "Events rejected, number of jets:" << nGoodJets;
0205     return false;
0206   }
0207
0208   double dEta = std::abs(genJetsWithoutLeptonsP4[0].eta() - genJetsWithoutLeptonsP4[1].eta());
0209   float invMassLeadingJet = (genJetsWithoutLeptonsP4[0] + genJetsWithoutLeptonsP4[1]).M();
0210
0211   if (dEta >= minDeltaEta && dEta <= maxDeltaEta && invMassLeadingJet > mininvmass) {
0212     return true;
0213   }
0214
0215   LogDebug("AJJJets") << "Events rejected, dEta:" << dEta << ", mjj:" << invMassLeadingJet;
0216   return false;
0217 }
0218
0219 void AJJGenJetFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0220   edm::ParameterSetDescription desc;
0221   desc.add<edm::InputTag>("GenJetCollection", edm::InputTag("ak4GenJetsNoNu"));
0222   desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
0223   desc.add<double>("minPt", -1.0)->setComment("If this is negative, no cut on jets is applied");
0224   desc.add<double>("minEta", -5.0);
0225   desc.add<double>("maxEta", 5.0);
0226   desc.add<double>("deltaRJetLep", 0.);
0227   desc.add<double>("minDeltaEta", 0.0);
0228   desc.add<double>("maxDeltaEta", 9999.0);
0229   desc.add<double>("MinInvMass", 0.0);
0230   desc.add<double>("maxPhotonEta", 5);
0231   desc.add<double>("minPhotonPt", 50);
0232   desc.add<double>("maxPhotonPt", 10000);
0233   descriptions.addDefault(desc);
0234 }
0235 //define this as a plug-in
0236 DEFINE_FWK_MODULE(AJJGenJetFilter);