AJJGenJet filter added.

GeneratorInterface/GenFilters/plugins/AJJGenJetFilter.cc

@@ -0,0 +1,157 @@
+#include "GeneratorInterface/GenFilters/plugins/AJJGenJetFilter.h"
+
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include "DataFormats/Math/interface/deltaR.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include <HepMC/GenVertex.h>
+
+// ROOT includes
+#include "TMath.h"
+#include "TH1.h"
+#include "TH2.h"
+
+// C++ includes
+#include <iostream>
+#include <vector>
+
+using namespace edm;
+using namespace std;
+
+AJJGenJetFilter::AJJGenJetFilter(const edm::ParameterSet& iConfig)
+    : ptMin(iConfig.getUntrackedParameter<double>("minPt", 0)),
+      etaMin(iConfig.getUntrackedParameter<double>("minEta", -4.5)),
+      etaMax(iConfig.getUntrackedParameter<double>("maxEta", 4.5)),
+      minDeltaEta(iConfig.getUntrackedParameter<double>("minDeltaEta", 0.0)),
+      maxDeltaEta(iConfig.getUntrackedParameter<double>("maxDeltaEta", 99999.0)),
+      deltaRJetLep(iConfig.getUntrackedParameter<double>("deltaRJetLep", 0.0)),
+      maxPhotonEta(iConfig.getUntrackedParameter<double>("maxPhotonEta", 5)),
+      minPhotonPt(iConfig.getUntrackedParameter<double>("minPhotonPt", 50)),
+      maxPhotonPt(iConfig.getUntrackedParameter<double>("maxPhotonPt", 10000)),
+      mininvmass(iConfig.getUntrackedParameter<double>("MinInvMass", 0.0)) {
+  m_GenJetCollection = consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("GenJetCollection"));
+  m_GenParticleCollection = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticles"));
+
+  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() {}
+
+vector<const reco::GenParticle*> AJJGenJetFilter::filterGenLeptons(const vector<reco::GenParticle>* particles) {
+  vector<const reco::GenParticle*> 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;
+}
+
+vector<const reco::GenParticle*> AJJGenJetFilter::filterGenPhotons(const vector<reco::GenParticle>* particles) {
+  vector<const reco::GenParticle*> out;
+
+  for (const auto& p : *particles) {
+    int absPdgId = std::abs(p.pdgId());
+
+    if ((absPdgId == 22) && p.isHardProcess()) {
+      if (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;
+}
+
+vector<const reco::GenJet*> AJJGenJetFilter::filterGenJets(const vector<reco::GenJet>* jets) {
+  vector<const reco::GenJet*> out;
+
+  for (unsigned i = 0; i < jets->size(); i++) {
+    const reco::GenJet* j = &((*jets)[i]);
+
+    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::Event& iEvent, const edm::EventSetup& iSetup) {
+  using namespace edm;
+
+  Handle<vector<reco::GenJet> > handleGenJets;
+  iEvent.getByToken(m_GenJetCollection, handleGenJets);
+  const vector<reco::GenJet>* genJets = handleGenJets.product();
+  // Getting filtered generator jets
+  vector<const reco::GenJet*> filGenJets = filterGenJets(genJets);
+
+  Handle<reco::GenParticleCollection> genParticelesCollection;
+  iEvent.getByToken(m_GenParticleCollection, genParticelesCollection);
+  const vector<reco::GenParticle>* genParticles = genParticelesCollection.product();
+  vector<const reco::GenParticle*> filGenLep = filterGenLeptons(genParticles);
+
+  // Getting p4 of jet with no lepton
+  vector<math::XYZTLorentzVector> genJetsWithoutLeptonsP4;
+  unsigned int jetIdx = 0;
+  unsigned int nGoodJets = 0;
+  while (jetIdx < filGenJets.size()) {
+    bool jetWhitoutLep = true;
+
+    const math::XYZTLorentzVector& p4J = (filGenJets[jetIdx])->p4();
+    for (unsigned int i = 0; i < filGenLep.size() && jetWhitoutLep; ++i) {
+      if (reco::deltaR2((filGenLep[i])->p4(), p4J) < deltaRJetLep * deltaRJetLep)
+        jetWhitoutLep = false;
+    }
+    if (jetWhitoutLep) {
+      if (genJetsWithoutLeptonsP4.size() < 2) {
+        genJetsWithoutLeptonsP4.push_back(p4J);
+      }
+      nGoodJets++;
+    };
+    ++jetIdx;
+  }
+
+  vector<const reco::GenParticle*> filGenPhotons = filterGenPhotons(genParticles);
+
+  if (filGenPhotons.size() != 1) {
+    edm::LogInfo("AJJPhoton") << "Events rejected, number of photons:" << filGenPhotons.size();
+    return false;
+  }
+
+  if (ptMin < 0)
+    return true;
+
+  //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 = fabs(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;
+}
+//define this as a plug-in
+DEFINE_FWK_MODULE(AJJGenJetFilter);

GeneratorInterface/GenFilters/plugins/AJJGenJetFilter.h

@@ -0,0 +1,68 @@
+#ifndef AJJGenJetAnalyzer_h
+#define AJJGenJetAnalyzer_h
+
+// CMSSW include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/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 "DataFormats/JetReco/interface/GenJetCollection.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
+
+// ROOT includes
+#include "TFile.h"
+#include "TH1D.h"
+#include "TH2D.h"
+
+// C++ include files
+#include <memory>
+#include <map>
+#include <vector>
+
+using namespace edm;
+using namespace std;
+//
+// class declaration
+//
+
+class AJJGenJetFilter : public edm::EDFilter {
+public:
+  explicit AJJGenJetFilter(const edm::ParameterSet& pset);
+  ~AJJGenJetFilter() override;
+
+  //void analyze(edm::Event&, const edm::EventSetup&)
+  bool filter(edm::Event&, const edm::EventSetup&) override;
+
+private:
+  // ----------memeber function----------------------
+  std::vector<const reco::GenJet*> filterGenJets(const vector<reco::GenJet>* jets);
+  std::vector<const reco::GenParticle*> filterGenLeptons(const std::vector<reco::GenParticle>* particles);
+  std::vector<const reco::GenParticle*> filterGenPhotons(const std::vector<reco::GenParticle>* particles);
+
+  //**************************
+  // Private Member data *****
+private:
+  // Dijet cut
+  double ptMin;
+  double etaMin;
+  double etaMax;
+  double minDeltaEta;
+  double maxDeltaEta;
+  double deltaRJetLep;
+  double maxPhotonEta;
+  double minPhotonPt;
+  double maxPhotonPt;
+  double mininvmass;
+
+  // Input tags
+  edm::EDGetTokenT<reco::GenJetCollection> m_GenJetCollection;
+  edm::EDGetTokenT<reco::GenParticleCollection> m_GenParticleCollection;
+};
+
+#endif

GeneratorInterface/GenFilters/python/AJJGenJetFilter_cff.py

@@ -0,0 +1,12 @@
+import FWCore.ParameterSet.Config as cms
+
+from GeneratorInterface.GenFilters.VJJGenJetFilter_cfi import *
+from RecoJets.Configuration.GenJetParticles_cff import *
+from RecoJets.Configuration.RecoGenJets_cff import *
+
+vjjGenJetFilterSeq = cms.Sequence(
+  genParticlesForJetsNoNu*
+  ak4GenJetsNoNu*
+  vjjGenJetFilterPhotonInBarrelMjj300
+)
+

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 )
+)