adding VBF MJJ enrichment for DY

cms-sw · Sep 5, 2017 · c1066c2 · c1066c2
1 parent 28b86df
commit c1066c2
Show file tree

Hide file tree

Showing 2 changed files with 69 additions and 2 deletions.
diff --git a/GeneratorInterface/GenFilters/interface/VBFGenJetFilter.h b/GeneratorInterface/GenFilters/interface/VBFGenJetFilter.h
@@ -12,6 +12,7 @@
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
 
 #include "DataFormats/JetReco/interface/GenJetCollection.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 
 // ROOT includes
 #include "TFile.h"
@@ -42,7 +43,7 @@ class VBFGenJetFilter : public edm::EDFilter {
   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);
+  std::vector<const reco::GenParticle*> filterGenLeptons(const std::vector<reco::GenParticle>* particles);
 
   //**************************
   // Private Member data *****
@@ -52,6 +53,7 @@ class VBFGenJetFilter : public edm::EDFilter {
 
   // Dijet cut
   bool   oppositeHemisphere;
+  bool   leadJetsNoLepMass;
   double ptMin;
   double etaMin;
   double etaMax;
@@ -61,9 +63,13 @@ class VBFGenJetFilter : public edm::EDFilter {
   double maxDeltaPhi;
   double minDeltaEta;
   double maxDeltaEta;
+  double minLeadingJetsInvMass;
+  double maxLeadingJetsInvMass;
+  double deltaRJetLep;
 
   // Input tags
   edm::EDGetTokenT< reco::GenJetCollection > m_inputTag_GenJetCollection;
+  edm::EDGetTokenT< reco::GenParticleCollection > m_inputTag_GenParticleCollection;
 
 
 };

diff --git a/GeneratorInterface/GenFilters/src/VBFGenJetFilter.cc b/GeneratorInterface/GenFilters/src/VBFGenJetFilter.cc
@@ -16,6 +16,7 @@ using namespace std;
 
 VBFGenJetFilter::VBFGenJetFilter(const edm::ParameterSet& iConfig) :
 oppositeHemisphere(iConfig.getUntrackedParameter<bool>  ("oppositeHemisphere",false)),
+leadJetsNoLepMass (iConfig.getUntrackedParameter<bool>  ("leadJetsNoLepMass", false)),
 ptMin             (iConfig.getUntrackedParameter<double>("minPt",                20)),
 etaMin            (iConfig.getUntrackedParameter<double>("minEta",             -5.0)),
 etaMax            (iConfig.getUntrackedParameter<double>("maxEta",              5.0)),
@@ -24,17 +25,41 @@ maxInvMass        (iConfig.getUntrackedParameter<double>("maxInvMass",      9999
 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))
+maxDeltaEta       (iConfig.getUntrackedParameter<double>("maxDeltaEta",     99999.0)),
+minLeadingJetsInvMass (iConfig.getUntrackedParameter<double>("minLeadingJetsInvMass",          0.0)),
+maxLeadingJetsInvMass (iConfig.getUntrackedParameter<double>("maxLeadingJetsInvMass",      99999.0)),
+deltaRJetLep          (iConfig.getUntrackedParameter<double>("deltaRJetLep",                   0.3))
 {
 
   m_inputTag_GenJetCollection = consumes<reco::GenJetCollection>(iConfig.getUntrackedParameter<edm::InputTag>("inputTag_GenJetCollection",edm::InputTag("ak5GenJetsNoNu")));
+  m_inputTag_GenParticleCollection  = consumes<reco::GenParticleCollection>(iConfig.getUntrackedParameter<edm::InputTag>("genParticles",edm::InputTag("genParticles")));
 
 }
 
 VBFGenJetFilter::~VBFGenJetFilter(){
 
 }
 
+
+vector<const reco::GenParticle*> VBFGenJetFilter::filterGenLeptons(const vector<reco::GenParticle>* particles){
+  vector<const reco::GenParticle*> out;
+
+  for(unsigned i=0; i<particles->size(); i++){
+
+      const reco::GenParticle* p = &((*particles)[i]);
+      int absPdgId = abs(p->pdgId());
+
+      if(((absPdgId == 11) || (absPdgId == 13) || (absPdgId == 15)) && p->isHardProcess()) {
+          out.push_back(p);              
+      }
+
+
+  }      
+  return out;
+}
+
+
+
 vector<const reco::GenJet*> VBFGenJetFilter::filterGenJets(const vector<reco::GenJet>* jets){
 
   vector<const reco::GenJet*> out;
@@ -68,6 +93,42 @@ bool VBFGenJetFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
   // If we do not find at least 2 jets veto the event
   if(filGenJets.size()<2){return false;}
 
+
+  // Testing dijet mass   
+  if(leadJetsNoLepMass) { 
+
+      Handle<reco::GenParticleCollection> genParticelesCollection;
+      iEvent.getByToken(m_inputTag_GenParticleCollection, genParticelesCollection);
+      const vector<reco::GenParticle>*  genParticles = genParticelesCollection.product();
+
+      // Getting filtered generator muons
+      vector<const reco::GenParticle*> filGenLep = filterGenLeptons(genParticles);
+
+      // Getting p4 of jet with no lepton
+      vector<math::XYZTLorentzVector> genJetsWithoutLeptonsP4;
+      unsigned int jetIdx = 0;
+
+      while(genJetsWithoutLeptonsP4.size()<2 && jetIdx < filGenJets.size()) {
+          bool jetWhitoutLep = true;
+          math::XYZTLorentzVector p4J= (filGenJets[jetIdx])->p4();
+          for(unsigned int i = 0; i < filGenLep.size() && jetWhitoutLep; ++i) {
+              float dPhi = reco::deltaPhi((filGenLep[i])->p4().phi(), p4J.phi());
+              float dEta = (filGenLep[i])->p4().eta()-p4J.eta();
+              if(dPhi*dPhi + dEta*dEta < deltaRJetLep*deltaRJetLep)
+                  jetWhitoutLep = false;
+          }
+          if (jetWhitoutLep)  genJetsWithoutLeptonsP4.push_back(p4J);
+          ++jetIdx;
+      }
+
+      // Checking the invariant mass of the leading jets
+      if (genJetsWithoutLeptonsP4.size() < 2) return false;
+      float invMassLeadingJet = (genJetsWithoutLeptonsP4[0] + genJetsWithoutLeptonsP4[1]).M();
+      if ( invMassLeadingJet > minLeadingJetsInvMass  && invMassLeadingJet < maxLeadingJetsInvMass) return true;
+      else return false;
+  }
+
+
   for(unsigned a=0; a<filGenJets.size(); a++){
     for(unsigned b=a+1; b<filGenJets.size(); b++){