VBF filters for QCD sample requested by VBF InvHiggs group

GeneratorInterface/GenFilters/interface/VBFGenJetFilter.h

@@ -0,0 +1,71 @@
+#ifndef VBFGenJetFilter_h
+#define VBFGenJetFilter_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"
+
+// ROOT includes
+#include "TFile.h"
+#include "TH1D.h"
+
+// C++ include files
+#include <memory>
+#include <map>
+
+//
+// class declaration
+//
+
+class VBFGenJetFilter : public edm::EDFilter {
+public:
+  explicit VBFGenJetFilter(const edm::ParameterSet&);
+  ~VBFGenJetFilter();
+
+  virtual bool filter(edm::Event&, const edm::EventSetup&);
+private:
+
+  // ----------memeber function----------------------
+  int    charge      (const int& Id);
+  std::vector<HepMC::GenParticle*> getVisibleDecayProducts(HepMC::GenParticle* particle);
+  std::vector<HepMC::GenParticle*> getNu (const HepMC::GenEvent* particles);
+  std::vector<HepMC::GenParticle*> getSt3(const HepMC::GenEvent* particles);
+  void printGenVector(std::vector<HepMC::GenParticle*> vec);
+  double nuMET(std::vector<HepMC::GenParticle*> vNu);
+
+  std::vector<const reco::GenJet*> filterGenJets(const std::vector<reco::GenJet>* jets);
+  //   std::vector<const reco::GenJet*> filterGenJets(const std::vector<reco::GenJet>* jets);
+
+  //**************************
+  // Private Member data *****
+private:
+
+
+
+  // Dijet cut
+  bool   oppositeHemisphere;
+  double ptMin;
+  double etaMin;
+  double etaMax;
+  double minInvMass;
+  double maxInvMass;
+  double minDeltaPhi;
+  double maxDeltaPhi;
+  double minDeltaEta;
+  double maxDeltaEta;
+
+  // Input tags
+  edm::EDGetTokenT< reco::GenJetCollection > m_inputTag_GenJetCollection;
+
+
+};
+
+#endif

GeneratorInterface/GenFilters/python/VBFGenJetFilter_cff.py

@@ -0,0 +1,23 @@
+import FWCore.ParameterSet.Config as cms
+
+from GeneratorInterface.GenFilters.VBFGenJetFilter_cfi import *
+from RecoJets.Configuration.GenJetParticles_cff import *
+from RecoJets.Configuration.RecoGenJets_cff import *
+
+vbfGenJetFilterASeq = cms.Sequence(
+  genParticlesForJetsNoNu*
+  ak4GenJetsNoNu*
+  vbfGenJetFilterA
+)
+
+vbfGenJetFilterBSeq = cms.Sequence(
+  genParticlesForJetsNoNu*
+  ak4GenJetsNoNu*
+  vbfGenJetFilterB
+)
+
+vbfGenJetFilterCSeq = cms.Sequence(
+  genParticlesForJetsNoNu*
+  ak4GenJetsNoNu*
+  vbfGenJetFilterC
+)

GeneratorInterface/GenFilters/python/VBFGenJetFilter_cfi.py

@@ -0,0 +1,52 @@
+import FWCore.ParameterSet.Config as cms
+
+vbfGenJetFilterA = cms.EDFilter("VBFGenJetFilter",
+
+  inputTag_GenJetCollection = cms.untracked.InputTag('ak4GenJetsNoNu'),
+
+  oppositeHemisphere = cms.untracked.bool  ( False), # Require j1_eta*j2_eta<0
+  minPt              = cms.untracked.double(    40), # Minimum dijet jet_pt
+  minEta             = cms.untracked.double(  -4.8), # Minimum dijet jet_eta
+  maxEta             = cms.untracked.double(   4.8), # Maximum dijet jet_eta
+  minInvMass         = cms.untracked.double( 1000.), # Minimum dijet invariant mass
+  maxInvMass         = cms.untracked.double(99999.), # Maximum dijet invariant mass
+  minDeltaPhi        = cms.untracked.double(  -1.0), # Minimum dijet delta phi
+  maxDeltaPhi        = cms.untracked.double(  2.15), # Maximum dijet delta phi
+  minDeltaEta        = cms.untracked.double(   3.0), # Minimum dijet delta eta
+  maxDeltaEta        = cms.untracked.double(99999.)  # Maximum dijet delta eta
+
+)
+
+vbfGenJetFilterB = cms.EDFilter("VBFGenJetFilter",
+
+  inputTag_GenJetCollection = cms.untracked.InputTag('ak4GenJetsNoNu'),
+
+  oppositeHemisphere = cms.untracked.bool  ( False), # Require j1_eta*j2_eta<0
+  minPt              = cms.untracked.double(    40), # Minimum dijet jet_pt
+  minEta             = cms.untracked.double(  -4.8), # Minimum dijet jet_eta
+  maxEta             = cms.untracked.double(   4.8), # Maximum dijet jet_eta
+  minInvMass         = cms.untracked.double( 1000.), # Minimum dijet invariant mass
+  maxInvMass         = cms.untracked.double(99999.), # Maximum dijet invariant mass
+  minDeltaPhi        = cms.untracked.double(  2.15), # Minimum dijet delta phi
+  maxDeltaPhi        = cms.untracked.double(   3.2), # Maximum dijet delta phi
+  minDeltaEta        = cms.untracked.double(   3.0), # Minimum dijet delta eta
+  maxDeltaEta        = cms.untracked.double(99999.)  # Maximum dijet delta eta
+
+)
+
+vbfGenJetFilterC = cms.EDFilter("VBFGenJetFilter",
+
+  inputTag_GenJetCollection = cms.untracked.InputTag('ak4GenJetsNoNu'),
+
+  oppositeHemisphere = cms.untracked.bool  ( False), # Require j1_eta*j2_eta<0
+  minPt              = cms.untracked.double(    40), # Minimum dijet jet_pt
+  minEta             = cms.untracked.double(  -4.8), # Minimum dijet jet_eta
+  maxEta             = cms.untracked.double(   4.8), # Maximum dijet jet_eta
+  minInvMass         = cms.untracked.double( 1000.), # Minimum dijet invariant mass
+  maxInvMass         = cms.untracked.double(99999.), # Maximum dijet invariant mass
+  minDeltaPhi        = cms.untracked.double(  -1.0), # Minimum dijet delta phi
+  maxDeltaPhi        = cms.untracked.double(   3.2), # Maximum dijet delta phi
+  minDeltaEta        = cms.untracked.double(   3.0), # Minimum dijet delta eta
+  maxDeltaEta        = cms.untracked.double(99999.)  # Maximum dijet delta eta
+
+)

GeneratorInterface/GenFilters/src/VBFGenJetFilter.cc

@@ -0,0 +1,104 @@
+#include "GeneratorInterface/GenFilters/interface/VBFGenJetFilter.h"
+
+#include "DataFormats/Math/interface/deltaPhi.h"
+
+#include <HepMC/GenVertex.h>
+
+// ROOT includes
+#include "TMath.h"
+
+// C++ includes
+#include <iostream>
+
+using namespace edm;
+using namespace std;
+
+
+VBFGenJetFilter::VBFGenJetFilter(const edm::ParameterSet& iConfig) :
+oppositeHemisphere(iConfig.getUntrackedParameter<bool>  ("oppositeHemisphere",false)),
+ptMin             (iConfig.getUntrackedParameter<double>("minPt",                20)),
+etaMin            (iConfig.getUntrackedParameter<double>("minEta",             -5.0)),
+etaMax            (iConfig.getUntrackedParameter<double>("maxEta",              5.0)),
+minInvMass        (iConfig.getUntrackedParameter<double>("minInvMass",          0.0)),
+maxInvMass        (iConfig.getUntrackedParameter<double>("maxInvMass",      99999.0)),
+minDeltaPhi       (iConfig.getUntrackedParameter<double>("minDeltaPhi",        -1.0)),
+maxDeltaPhi       (iConfig.getUntrackedParameter<double>("maxDeltaPhi",     99999.0)),
+minDeltaEta       (iConfig.getUntrackedParameter<double>("minDeltaEta",        -1.0)),
+maxDeltaEta       (iConfig.getUntrackedParameter<double>("maxDeltaEta",     99999.0))
+{
+
+  m_inputTag_GenJetCollection = consumes<reco::GenJetCollection>(iConfig.getUntrackedParameter<edm::InputTag>("inputTag_GenJetCollection",edm::InputTag("ak5GenJetsNoNu")));
+
+}
+
+VBFGenJetFilter::~VBFGenJetFilter(){
+
+}
+
+vector<const reco::GenJet*> VBFGenJetFilter::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);
+    }
+  }
+
+  return out;
+}
+
+
+// ------------ method called to skim the data  ------------
+bool VBFGenJetFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
+{
+  using namespace edm;
+
+  Handle< vector<reco::GenJet> > handleGenJets;
+  iEvent.getByToken(m_inputTag_GenJetCollection, handleGenJets);
+  const vector<reco::GenJet>* genJets = handleGenJets.product();
+
+  // Getting filtered generator jets
+  vector<const reco::GenJet*> filGenJets = filterGenJets(genJets);
+
+  // If we do not find at least 2 jets veto the event
+  if(filGenJets.size()<2){return false;}
+
+  for(unsigned a=0; a<filGenJets.size(); a++){
+    for(unsigned b=a+1; b<filGenJets.size(); b++){    
+
+      const reco::GenJet* pA = filGenJets[a];
+      const reco::GenJet* pB = filGenJets[b];
+
+      // Getting the dijet vector
+      math::XYZTLorentzVector diJet = pA->p4() + pB->p4();
+
+      // Testing opposite hemispheres
+      double dijetProd = pA->p4().eta()*pB->p4().eta();
+      if(oppositeHemisphere && dijetProd>=0){continue;}
+
+      // Testing dijet mass
+      double invMass = diJet.mass();
+      if(invMass<=minInvMass || invMass>maxInvMass){continue;}
+
+      // Testing dijet delta eta
+      double dEta = fabs(pA->p4().eta()-pB->p4().eta());
+      if(dEta<=minDeltaEta || dEta>maxDeltaEta){continue;}
+
+      // Testing dijet delta phi
+      double dPhi = fabs(reco::deltaPhi(pA->p4().phi(),pB->p4().phi()));
+      if(dPhi<=minDeltaPhi || dPhi>maxDeltaPhi){continue;}
+
+      return true;
+    }
+  }
+
+  return false;
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(VBFGenJetFilter);