Skip to content

Commit

Permalink
Merge pull request #20392 from gmandorl/modified_VBF_filter
Browse files Browse the repository at this point in the history 
Modified vbf filter
  • Loading branch information
@cmsbuild
cmsbuild committed Sep 19, 2017
2 parents 4b2f12a + 9b919ff commit bd6cf02
Show file tree
Hide file tree
Showing 4 changed files with 202 additions and 15 deletions.
1 change: 1 addition & 0 deletions GeneratorInterface/GenFilters/BuildFile.xml
View file
Expand Up @@ -25,4 +25,5 @@
<use name="DataFormats/HepMCCandidate"/>
<use name="DataFormats/JetReco"/>
<use name="DataFormats/EgammaReco"/>
<use name="DataFormats/Math"/>
<flags EDM_PLUGIN="1"/>
19 changes: 17 additions & 2 deletions GeneratorInterface/GenFilters/interface/VBFGenJetFilter.h
View file
Expand Up @@ -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"
Expand All @@ -21,6 +22,7 @@
#include <memory>
#include <map>

#include <unordered_set>
//
// class declaration
//
Expand All @@ -41,8 +43,13 @@ class VBFGenJetFilter : public edm::EDFilter {
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);
std::vector<const reco::GenJet*> filterGenJets(const std::vector<reco::GenJet>* jets);
std::vector<const reco::GenParticle*> filterGenLeptons(const std::vector<reco::GenParticle>* particles);

bool isHardProcess(const reco::GenParticle &p);
const reco::GenParticle * firstCopy(const reco::GenParticle &p);
const reco::GenParticle * previousCopy(const reco::GenParticle &p);
const reco::GenParticle * nextCopy(const reco::GenParticle &p);

//**************************
// Private Member data *****
Expand All @@ -52,20 +59,28 @@ class VBFGenJetFilter : public edm::EDFilter {

// Dijet cut
bool oppositeHemisphere;
bool leadJetsNoLepMass;
double ptMin;
double etaMin;
double etaMax;
double minInvMass;
double maxInvMass;
double minLeadingJetsInvMass;
double maxLeadingJetsInvMass;
double deltaRNoLep;
double minDeltaPhi;
double maxDeltaPhi;
double minDeltaEta;
double maxDeltaEta;

// Input tags
edm::EDGetTokenT< reco::GenJetCollection > m_inputTag_GenJetCollection;
edm::EDGetTokenT< reco::GenParticleCollection > m_inputTag_GenParticleCollection;


};




#endif
14 changes: 13 additions & 1 deletion GeneratorInterface/GenFilters/python/VBFGenJetFilter_cfi.py
View file
Expand Up @@ -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), # Maximum dijet invariant mass

)
183 changes: 171 additions & 12 deletions GeneratorInterface/GenFilters/src/VBFGenJetFilter.cc
View file
@@ -1,5 +1,6 @@
#include "GeneratorInterface/GenFilters/interface/VBFGenJetFilter.h"

#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include <HepMC/GenVertex.h>
Expand All @@ -15,26 +16,132 @@ 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))
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)),
minInvMass (iConfig.getUntrackedParameter<double>("minInvMass", 0.0)),
maxInvMass (iConfig.getUntrackedParameter<double>("maxInvMass", 99999.0)),
minLeadingJetsInvMass (iConfig.getUntrackedParameter<double>("minLeadingJetsInvMass", 0.0)),
maxLeadingJetsInvMass (iConfig.getUntrackedParameter<double>("maxLeadingJetsInvMass", 99999.0)),
deltaRNoLep (iConfig.getUntrackedParameter<double>("deltaRNoLep", 0.3)),
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")));

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")));
// m_inputTag_GenParticleCollection = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticles"));

}

VBFGenJetFilter::~VBFGenJetFilter(){

}


bool VBFGenJetFilter::isHardProcess(const reco::GenParticle &p) {

//status 3 in pythia6 means hard process;
if (p.status()==3) return true;

//hard process codes for pythia8 are 21-29 inclusive (currently 21,22,23,24 are used)
if (p.status()>20 && p.status()<30) return true;

//if this is a final state or decayed particle,
//check if direct mother is a resonance decay in pythia8 but exclude FSR branchings
//(In pythia8 if a resonance decay product did not undergo any further branchings
//it will be directly stored as status 1 or 2 without any status 23 copy)
if (p.status()==1 || p.status()==2) {
// const reco::GenParticle *um = mother(p);
const reco::GenParticle *um = static_cast<const reco::GenParticle*>(p.mother(0));
if (um) {
const reco::GenParticle *firstcopy = firstCopy(*um);
bool fromResonance = firstcopy && firstcopy->status()==22;

const reco::GenParticle *umNext = nextCopy(*um);
bool fsrBranching = umNext && umNext->status()>50 && umNext->status()<60;

if (fromResonance && !fsrBranching) return true;
}
}

return false;

}

/////////////////////////////////////////////////////////////////////////////

const reco::GenParticle * VBFGenJetFilter::firstCopy(const reco::GenParticle &p) {
const reco::GenParticle *pcopy = &p;
std::unordered_set<const reco::GenParticle*> dupCheck;
while (previousCopy(*pcopy)) {
dupCheck.insert(pcopy);
pcopy = previousCopy(*pcopy);
if (dupCheck.count(pcopy)) return 0;
}
return pcopy;
}

/////////////////////////////////////////////////////////////////////////////

const reco::GenParticle * VBFGenJetFilter::previousCopy(const reco::GenParticle &p) {

const unsigned int nmoth = p.numberOfMothers();
for (unsigned int imoth = 0; imoth<nmoth; ++imoth) {
const reco::GenParticle *moth = static_cast<const reco::GenParticle*>(p.mother(imoth));//mother(p,imoth);
if (moth->pdgId()==p.pdgId()) {
return moth;
}
}

return 0;
}

/////////////////////////////////////////////////////////////////////////////

const reco::GenParticle * VBFGenJetFilter::nextCopy(const reco::GenParticle &p) {

const unsigned int ndau = p.numberOfDaughters();
for (unsigned int idau = 0; idau<ndau; ++idau) {
const reco::GenParticle *dau = static_cast<const reco::GenParticle*>(p.daughter(idau));//daughter(p,idau);
if (dau->pdgId()==p.pdgId()) {
return dau;
}
}

return 0;
}

/////////////////////////////////////////////////////////////////////////////

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)) && isHardProcess(p)) {
out.push_back(&p);
}


}
return out;
}


/////////////////////////////////////////////////////////////////////////////


vector<const reco::GenJet*> VBFGenJetFilter::filterGenJets(const vector<reco::GenJet>* jets){

vector<const reco::GenJet*> out;
Expand All @@ -53,6 +160,11 @@ vector<const reco::GenJet*> VBFGenJetFilter::filterGenJets(const vector<reco::Ge
}


/////////////////////////////////////////////////////////////////////////////




// ------------ method called to skim the data ------------
bool VBFGenJetFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
Expand All @@ -62,12 +174,58 @@ bool VBFGenJetFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
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;}


// 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) < deltaRNoLep*deltaRNoLep)

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

Expand All @@ -85,6 +243,7 @@ bool VBFGenJetFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
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;}
Expand Down

0 comments on commit bd6cf02

Please sign in to comment.