GenFilters/plugins/VBFGenJetFilter.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0003 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 #include "DataFormats/Math/interface/deltaR.h"
0006 #include "DataFormats/Math/interface/LorentzVector.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/global/EDFilter.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/EDGetToken.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014
0015 #include <cmath>
0016 #include <cstdlib>
0017 #include <vector>
0018
0019 class VBFGenJetFilter : public edm::global::EDFilter<> {
0020 public:
0021   explicit VBFGenJetFilter(const edm::ParameterSet&);
0022
0023   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0024
0025 private:
0026   std::vector<const reco::GenJet*> filterGenJets(const std::vector<reco::GenJet>* jets) const;
0027   std::vector<const reco::GenParticle*> filterGenLeptons(const std::vector<reco::GenParticle>* particles) const;
0028
0029   // Dijet cut
0030   const bool oppositeHemisphere;
0031   const bool leadJetsNoLepMass;
0032   const double ptMin;
0033   const double etaMin;
0034   const double etaMax;
0035   const double minInvMass;
0036   const double maxInvMass;
0037   const double minDeltaPhi;
0038   const double maxDeltaPhi;
0039   const double minDeltaEta;
0040   const double maxDeltaEta;
0041   const double minLeadingJetsInvMass;
0042   const double maxLeadingJetsInvMass;
0043   const double deltaRJetLep;
0044
0045   // Input tags
0046   edm::EDGetTokenT<reco::GenJetCollection> m_inputTag_GenJetCollection;
0047   edm::EDGetTokenT<reco::GenParticleCollection> m_inputTag_GenParticleCollection;
0048 };
0049
0050 using namespace std;
0051
0052 VBFGenJetFilter::VBFGenJetFilter(const edm::ParameterSet& iConfig)
0053     : oppositeHemisphere(iConfig.getUntrackedParameter<bool>("oppositeHemisphere", false)),
0054       leadJetsNoLepMass(iConfig.getUntrackedParameter<bool>("leadJetsNoLepMass", false)),
0055       ptMin(iConfig.getUntrackedParameter<double>("minPt", 20)),
0056       etaMin(iConfig.getUntrackedParameter<double>("minEta", -5.0)),
0057       etaMax(iConfig.getUntrackedParameter<double>("maxEta", 5.0)),
0058       minInvMass(iConfig.getUntrackedParameter<double>("minInvMass", 0.0)),
0059       maxInvMass(iConfig.getUntrackedParameter<double>("maxInvMass", 99999.0)),
0060       minDeltaPhi(iConfig.getUntrackedParameter<double>("minDeltaPhi", -1.0)),
0061       maxDeltaPhi(iConfig.getUntrackedParameter<double>("maxDeltaPhi", 99999.0)),
0062       minDeltaEta(iConfig.getUntrackedParameter<double>("minDeltaEta", -1.0)),
0063       maxDeltaEta(iConfig.getUntrackedParameter<double>("maxDeltaEta", 99999.0)),
0064       minLeadingJetsInvMass(iConfig.getUntrackedParameter<double>("minLeadingJetsInvMass", 0.0)),
0065       maxLeadingJetsInvMass(iConfig.getUntrackedParameter<double>("maxLeadingJetsInvMass", 99999.0)),
0066       deltaRJetLep(iConfig.getUntrackedParameter<double>("deltaRJetLep", 0.3)) {
0067   m_inputTag_GenJetCollection = consumes<reco::GenJetCollection>(
0068       iConfig.getUntrackedParameter<edm::InputTag>("inputTag_GenJetCollection", edm::InputTag("ak5GenJetsNoNu")));
0069   if (leadJetsNoLepMass)
0070     m_inputTag_GenParticleCollection = consumes<reco::GenParticleCollection>(
0071         iConfig.getUntrackedParameter<edm::InputTag>("genParticles", edm::InputTag("genParticles")));
0072 }
0073
0074 vector<const reco::GenParticle*> VBFGenJetFilter::filterGenLeptons(const vector<reco::GenParticle>* particles) const {
0075   vector<const reco::GenParticle*> out;
0076
0077   for (const auto& p : *particles) {
0078     int absPdgId = std::abs(p.pdgId());
0079
0080     if (((absPdgId == 11) || (absPdgId == 13) || (absPdgId == 15)) && p.isHardProcess()) {
0081       out.push_back(&p);
0082     }
0083   }
0084   return out;
0085 }
0086
0087 vector<const reco::GenJet*> VBFGenJetFilter::filterGenJets(const vector<reco::GenJet>* jets) const {
0088   vector<const reco::GenJet*> out;
0089
0090   for (unsigned i = 0; i < jets->size(); i++) {
0091     const reco::GenJet* j = &((*jets)[i]);
0092
0093     if (j->p4().pt() > ptMin && j->p4().eta() > etaMin && j->p4().eta() < etaMax) {
0094       out.push_back(j);
0095     }
0096   }
0097
0098   return out;
0099 }
0100
0101 // ------------ method called to skim the data  ------------
0102 bool VBFGenJetFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0103   edm::Handle<vector<reco::GenJet> > handleGenJets;
0104   iEvent.getByToken(m_inputTag_GenJetCollection, handleGenJets);
0105   const vector<reco::GenJet>* genJets = handleGenJets.product();
0106
0107   // Getting filtered generator jets
0108   vector<const reco::GenJet*> filGenJets = filterGenJets(genJets);
0109
0110   // If we do not find at least 2 jets veto the event
0111   if (filGenJets.size() < 2) {
0112     return false;
0113   }
0114
0115   // Testing dijet mass
0116   if (leadJetsNoLepMass) {
0117     edm::Handle<reco::GenParticleCollection> genParticelesCollection;
0118     iEvent.getByToken(m_inputTag_GenParticleCollection, genParticelesCollection);
0119     const vector<reco::GenParticle>* genParticles = genParticelesCollection.product();
0120
0121     // Getting filtered generator muons
0122     vector<const reco::GenParticle*> filGenLep = filterGenLeptons(genParticles);
0123
0124     // Getting p4 of jet with no lepton
0125     vector<math::XYZTLorentzVector> genJetsWithoutLeptonsP4;
0126     unsigned int jetIdx = 0;
0127
0128     while (genJetsWithoutLeptonsP4.size() < 2 && jetIdx < filGenJets.size()) {
0129       bool jetWhitoutLep = true;
0130
0131       const math::XYZTLorentzVector& p4J = (filGenJets[jetIdx])->p4();
0132       for (unsigned int i = 0; i < filGenLep.size() && jetWhitoutLep; ++i) {
0133         if (reco::deltaR2((filGenLep[i])->p4(), p4J) < deltaRJetLep * deltaRJetLep)
0134           jetWhitoutLep = false;
0135       }
0136       if (jetWhitoutLep)
0137         genJetsWithoutLeptonsP4.push_back(p4J);
0138       ++jetIdx;
0139     }
0140
0141     // Checking the invariant mass of the leading jets
0142     if (genJetsWithoutLeptonsP4.size() < 2)
0143       return false;
0144     float invMassLeadingJet = (genJetsWithoutLeptonsP4[0] + genJetsWithoutLeptonsP4[1]).M();
0145     if (invMassLeadingJet > minLeadingJetsInvMass && invMassLeadingJet < maxLeadingJetsInvMass)
0146       return true;
0147     else
0148       return false;
0149   }
0150
0151   for (unsigned a = 0; a < filGenJets.size(); a++) {
0152     for (unsigned b = a + 1; b < filGenJets.size(); b++) {
0153       const reco::GenJet* pA = filGenJets[a];
0154       const reco::GenJet* pB = filGenJets[b];
0155
0156       // Getting the dijet vector
0157       math::XYZTLorentzVector diJet = pA->p4() + pB->p4();
0158
0159       // Testing opposite hemispheres
0160       double dijetProd = pA->p4().eta() * pB->p4().eta();
0161       if (oppositeHemisphere && dijetProd >= 0) {
0162         continue;
0163       }
0164
0165       // Testing dijet mass
0166       double invMass = diJet.mass();
0167       if (invMass <= minInvMass || invMass > maxInvMass) {
0168         continue;
0169       }
0170
0171       // Testing dijet delta eta
0172       double dEta = fabs(pA->p4().eta() - pB->p4().eta());
0173       if (dEta <= minDeltaEta || dEta > maxDeltaEta) {
0174         continue;
0175       }
0176
0177       // Testing dijet delta phi
0178       double dPhi = fabs(reco::deltaPhi(pA->p4().phi(), pB->p4().phi()));
0179       if (dPhi <= minDeltaPhi || dPhi > maxDeltaPhi) {
0180         continue;
0181       }
0182
0183       return true;
0184     }
0185   }
0186
0187   return false;
0188 }
0189
0190 //define this as a plug-in
0191 DEFINE_FWK_MODULE(VBFGenJetFilter);