File indexing completed on 2024-04-06 12:13:39
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
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
0042 #include <memory>
0043 #include <map>
0044 #include <vector>
0045 #include <iostream>
0046
0047 using namespace edm;
0048 using namespace std;
0049
0050
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
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
0069 private:
0070
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
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
0164 bool AJJGenJetFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0165 using namespace edm;
0166
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
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
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
0236 DEFINE_FWK_MODULE(AJJGenJetFilter);