Merge pull request #24651 from gmandorl/from-CMSSW_9_3_9_patch1

VBFGenJetFilter has been modfied in order to select generator level Mjj

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"
@@ -28,9 +29,9 @@
 class VBFGenJetFilter : public edm::EDFilter {
 public:
   explicit VBFGenJetFilter(const edm::ParameterSet&);
-  ~VBFGenJetFilter();
+  ~VBFGenJetFilter() override;
 
-  virtual bool filter(edm::Event&, const edm::EventSetup&);
+  bool filter(edm::Event&, const edm::EventSetup&) override;
 private:
 
   // ----------memeber function----------------------
@@ -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;
 
 
 };

GeneratorInterface/GenFilters/python/VBFGenJetFilter_cfi.py

@@ -49,4 +49,16 @@
   minDeltaEta        = cms.untracked.double(   3.0), # Minimum dijet delta eta
   maxDeltaEta        = cms.untracked.double(99999.)  # Maximum dijet delta eta
 
-)
+)
+
+
+vbfGenJetFilterD = cms.EDFilter("VBFGenJetFilter",
+
+  inputTag_GenJetCollection = cms.untracked.InputTag('ak4GenJetsNoNu'),
+
+  leadJetsNoLepMass         = cms.untracked.bool  ( True),  # Require the cut on the mass of the leading jets
+  minLeadingJetsInvMass     = cms.untracked.double(   0.0), # Minimum dijet invariant mass
+  maxLeadingJetsInvMass     = cms.untracked.double(99999.), # Maximum dijet invariant mass
+  deltaRNoLep               = cms.untracked.double(   0.3), # Minimum deltaR(lepton, jet)
+
+)

GeneratorInterface/GenFilters/src/VBFGenJetFilter.cc

@@ -1,6 +1,7 @@
 #include "GeneratorInterface/GenFilters/interface/VBFGenJetFilter.h"
 
 #include "DataFormats/Math/interface/deltaPhi.h"
+#include "DataFormats/Math/interface/deltaR.h"
 
 #include <HepMC/GenVertex.h>
 
@@ -16,6 +17,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 +26,42 @@ 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")));
+  if (leadJetsNoLepMass) 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(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::GenJet*> VBFGenJetFilter::filterGenJets(const vector<reco::GenJet>* jets){
 
   vector<const reco::GenJet*> out;
@@ -68,6 +95,41 @@ 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;
+
+          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)  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++){