Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-26 10:51:19

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   //Handle<reco::GenParticleCollection> genParticelesCollection;
0169   //iEvent.getByToken(m_GenParticleCollection, genParticelesCollection);
0170   //const vector<reco::GenParticle>* genParticles = genParticelesCollection.product();
0171   auto const& genParticles = iEvent.get(m_GenParticleCollection);
0172   vector<const reco::GenParticle*> filGenLep = filterGenLeptons(&genParticles);
0173 
0174   vector<const reco::GenParticle*> filGenPhotons = filterGenPhotons(&genParticles);
0175 
0176   if (filGenPhotons.size() != 1) {
0177     edm::LogInfo("AJJPhoton") << "Events rejected, number of photons:" << filGenPhotons.size();
0178     return false;
0179   }
0180 
0181   if (ptMin < 0)
0182     return true;
0183 
0184   //Handle<vector<reco::GenJet> > handleGenJets;
0185   //iEvent.getByToken(m_GenJetCollection, handleGenJets);
0186   //const vector<reco::GenJet>* genJets = handleGenJets.product();
0187   auto const& genJets = iEvent.get(m_GenJetCollection);
0188 
0189   vector<const reco::GenJet*> filGenJets = filterGenJets(&genJets);
0190   // Getting p4 of jet with no lepton
0191   vector<math::XYZTLorentzVector> genJetsWithoutLeptonsP4;
0192   unsigned int nGoodJets = 0;
0193   for (auto const& j : filGenJets) {
0194     bool cleanJet = true;
0195     const math::XYZTLorentzVector& p4J = j->p4();
0196     for (auto const& p : filGenLep)
0197       if (reco::deltaR2(p->p4(), p4J) < deltaRJetLep * deltaRJetLep)
0198         cleanJet = false;
0199 
0200     if (cleanJet) {
0201       if (genJetsWithoutLeptonsP4.size() < 2)
0202         genJetsWithoutLeptonsP4.push_back(p4J);
0203       nGoodJets++;
0204     }
0205   }
0206 
0207   //If we do not find at least 2 jets veto the event
0208   if (nGoodJets < 2) {
0209     edm::LogInfo("AJJJets") << "Events rejected, number of jets:" << nGoodJets;
0210     return false;
0211   }
0212 
0213   double dEta = std::abs(genJetsWithoutLeptonsP4[0].eta() - genJetsWithoutLeptonsP4[1].eta());
0214   float invMassLeadingJet = (genJetsWithoutLeptonsP4[0] + genJetsWithoutLeptonsP4[1]).M();
0215 
0216   if (dEta >= minDeltaEta && dEta <= maxDeltaEta && invMassLeadingJet > mininvmass) {
0217     return true;
0218   }
0219 
0220   edm::LogInfo("AJJJets") << "Events rejected, dEta:" << dEta << ", mjj:" << invMassLeadingJet;
0221   return false;
0222 }
0223 
0224 void AJJGenJetFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0225   edm::ParameterSetDescription desc;
0226   desc.add<edm::InputTag>("GenJetCollection", edm::InputTag("ak4GenJetsNoNu"));
0227   desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
0228   desc.add<double>("minPt", -1.0)->setComment("If this is negative, no cut on jets is applied");
0229   desc.add<double>("minEta", -5.0);
0230   desc.add<double>("maxEta", 5.0);
0231   desc.add<double>("deltaRJetLep", 0.);
0232   desc.add<double>("minDeltaEta", 0.0);
0233   desc.add<double>("maxDeltaEta", 9999.0);
0234   desc.add<double>("MinInvMass", 0.0);
0235   desc.add<double>("maxPhotonEta", 5);
0236   desc.add<double>("minPhotonPt", 50);
0237   desc.add<double>("maxPhotonPt", 10000);
0238   descriptions.addDefault(desc);
0239 }
0240 //define this as a plug-in
0241 DEFINE_FWK_MODULE(AJJGenJetFilter);