Skip to content

Commit f8a2a28

Browse files
ferencekferencek
committed
new jet flavour definition
1 parent 7d19af0 commit f8a2a28

18 files changed

+1598
-37
lines changed

‎PhysicsTools/JetExamples/test/BuildFile.xml

+7
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,12 @@
11
<use name="clhep"/>
22
<use name="DataFormats/BTauReco"/>
3+
<library file="printJetFlavourInfo.cc">
4+
<flags EDM_PLUGIN="1"/>
5+
<use name="FWCore/Framework"/>
6+
<use name="SimDataFormats/JetMatching"/>
7+
<use name="SimGeneral/HepPDTRecord"/>
8+
<use name="DataFormats/JetReco"/>
9+
</library>
310
<library file="printJetFlavour.cc">
411
<flags EDM_PLUGIN="1"/>
512
<use name="FWCore/Framework"/>

‎PhysicsTools/JetExamples/test/printJetFlavour.py

+13-9
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,11 @@
11
import FWCore.ParameterSet.Config as cms
22

33
process = cms.Process("testJET")
4+
5+
process.options = cms.untracked.PSet(
6+
wantSummary = cms.untracked.bool(True)
7+
)
8+
49
process.load("FWCore.MessageService.MessageLogger_cfi")
510

611
process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi")
@@ -9,25 +14,26 @@
914
input = cms.untracked.int32(5)
1015
)
1116
process.source = cms.Source("PoolSource",
12-
fileNames = cms.untracked.vstring('/store/relval/2008/5/20/RelVal-RelValTTbar-1211209682-FakeConditions-2nd/0000/08765709-5826-DD11-9CE8-000423D94700.root')
17+
fileNames = cms.untracked.vstring('/store/relval/CMSSW_5_3_11_patch3/RelValTTbar/GEN-SIM-RECO/PU_START53_LV2_09Jul2013-v1/00000/A26AB2F9-41E9-E211-806F-002590593920.root')
1318
)
1419

15-
process.printTree = cms.EDFilter("ParticleListDrawer",
20+
process.printTree = cms.EDAnalyzer("ParticleListDrawer",
1621
src = cms.InputTag("genParticles"),
1722
maxEventsToPrint = cms.untracked.int32(1)
1823
)
1924

20-
process.myPartons = cms.EDFilter("PartonSelector",
25+
process.myPartons = cms.EDProducer("PartonSelector",
26+
src = cms.InputTag("genParticles"),
2127
withLeptons = cms.bool(False)
2228
)
2329

24-
process.flavourByRef = cms.EDFilter("JetPartonMatcher",
25-
jets = cms.InputTag("iterativeCone5CaloJets"),
30+
process.flavourByRef = cms.EDProducer("JetPartonMatcher",
31+
jets = cms.InputTag("ak5PFJets"),
2632
coneSizeToAssociate = cms.double(0.3),
2733
partons = cms.InputTag("myPartons")
2834
)
2935

30-
process.flavourByVal = cms.EDFilter("JetFlavourIdentifier",
36+
process.flavourByVal = cms.EDProducer("JetFlavourIdentifier",
3137
srcByReference = cms.InputTag("flavourByRef"),
3238
physicsDefinition = cms.bool(False)
3339
)
@@ -38,10 +44,8 @@
3844
srcByValue = cms.InputTag("flavourByVal")
3945
)
4046

41-
process.printEventNumber = cms.OutputModule("AsciiOutputModule")
42-
4347
process.p = cms.Path(process.printTree*process.myPartons*process.flavourByRef*process.flavourByVal*process.printEvent)
44-
process.outpath = cms.EndPath(process.printEventNumber)
48+
4549
process.MessageLogger.destinations = cms.untracked.vstring('cout','cerr')
4650
#process.MessageLogger.cout = cms.PSet(
4751
# threshold = cms.untracked.string('ERROR')

‎PhysicsTools/JetExamples/test/printJetFlavourInfo.cc

+248
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,248 @@
1+
// user include files
2+
#include "FWCore/Framework/interface/Frameworkfwd.h"
3+
#include "FWCore/Framework/interface/EDAnalyzer.h"
4+
5+
#include "FWCore/Framework/interface/Event.h"
6+
#include "FWCore/Framework/interface/MakerMacros.h"
7+
8+
#include "FWCore/ParameterSet/interface/ParameterSet.h"
9+
10+
#include "DataFormats/JetReco/interface/Jet.h"
11+
#include "SimDataFormats/JetMatching/interface/JetFlavourInfo.h"
12+
#include "SimDataFormats/JetMatching/interface/JetFlavourInfoMatching.h"
13+
14+
#include "DataFormats/Math/interface/deltaR.h"
15+
16+
// system include files
17+
#include <memory>
18+
19+
20+
class printJetFlavourInfo : public edm::EDAnalyzer {
21+
public:
22+
explicit printJetFlavourInfo(const edm::ParameterSet & );
23+
~printJetFlavourInfo() {};
24+
void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
25+
26+
private:
27+
edm::InputTag jetFlavourInfos_;
28+
edm::InputTag subjetFlavourInfos_;
29+
edm::InputTag groomedJets_;
30+
bool useSubjets_;
31+
};
32+
33+
printJetFlavourInfo::printJetFlavourInfo(const edm::ParameterSet& iConfig)
34+
{
35+
jetFlavourInfos_ = iConfig.getParameter<edm::InputTag>("jetFlavourInfos");
36+
subjetFlavourInfos_ = ( iConfig.exists("subjetFlavourInfos") ? iConfig.getParameter<edm::InputTag>("subjetFlavourInfos") : edm::InputTag() );
37+
groomedJets_ = ( iConfig.exists("groomedJets") ? iConfig.getParameter<edm::InputTag>("groomedJets") : edm::InputTag() );
38+
useSubjets_ = ( iConfig.exists("subjetFlavourInfos") && iConfig.exists("groomedJets") );
39+
}
40+
41+
void printJetFlavourInfo::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
42+
{
43+
edm::Handle<reco::JetFlavourInfoMatchingCollection> theJetFlavourInfos;
44+
iEvent.getByLabel(jetFlavourInfos_, theJetFlavourInfos );
45+
46+
edm::Handle<reco::JetFlavourInfoMatchingCollection> theSubjetFlavourInfos;
47+
edm::Handle<edm::View<reco::Jet> > groomedJets;
48+
49+
std::vector<int> matchedIndices;
50+
if( useSubjets_ )
51+
{
52+
iEvent.getByLabel(subjetFlavourInfos_, theSubjetFlavourInfos);
53+
iEvent.getByLabel(groomedJets_, groomedJets);
54+
55+
// match groomed and original jet
56+
std::vector<bool> jetLocks(theJetFlavourInfos->size(),false);
57+
std::vector<int> jetIndices;
58+
59+
for(size_t gj=0; gj<groomedJets->size(); ++gj)
60+
{
61+
double matchedDR = 1e9;
62+
int matchedIdx = -1;
63+
64+
for(reco::JetFlavourInfoMatchingCollection::const_iterator j = theJetFlavourInfos->begin();
65+
j != theJetFlavourInfos->end();
66+
++j)
67+
{
68+
if( jetLocks.at(j - theJetFlavourInfos->begin()) ) continue; // skip jets that have already been matched
69+
70+
double tempDR = reco::deltaR( j->first->rapidity(), j->first->phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi() );
71+
if( tempDR < matchedDR )
72+
{
73+
matchedDR = tempDR;
74+
matchedIdx = (j - theJetFlavourInfos->begin());
75+
}
76+
}
77+
78+
if( matchedIdx>=0 ) jetLocks.at(matchedIdx) = true;
79+
jetIndices.push_back(matchedIdx);
80+
}
81+
82+
if( std::find( jetIndices.begin(), jetIndices.end(), -1 ) != jetIndices.end() )
83+
throw cms::Exception("Jet matching failed") << "Matching groomed to original jets failed. Please check that the jet algorithm, jet size, and Pt threshold match for the two jet collections.";
84+
85+
for(size_t j=0; j<theJetFlavourInfos->size(); ++j)
86+
{
87+
std::vector<int>::iterator matchedIndex = std::find( jetIndices.begin(), jetIndices.end(), j );
88+
89+
matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
90+
}
91+
}
92+
93+
for ( reco::JetFlavourInfoMatchingCollection::const_iterator j = theJetFlavourInfos->begin();
94+
j != theJetFlavourInfos->end();
95+
++j ) {
96+
std::cout << "-------------------- Jet Flavour Info --------------------" << std::endl;
97+
98+
const reco::Jet *aJet = (*j).first.get();
99+
reco::JetFlavourInfo aInfo = (*j).second;
100+
std::cout << std::setprecision(2) << std::setw(6) << std::fixed
101+
<< "[printJetFlavourInfo] Jet " << (j - theJetFlavourInfos->begin()) << " pt, eta, rapidity, phi = " << aJet->pt() << ", "
102+
<< aJet->eta() << ", "
103+
<< aJet->rapidity() << ", "
104+
<< aJet->phi() << ", "
105+
<< std::endl;
106+
// ----------------------- Hadrons -------------------------------
107+
std::cout << " Hadron-based flavour: " << aInfo.getHadronFlavour() << std::endl;
108+
109+
const reco::GenParticleRefVector & bHadrons = aInfo.getbHadrons();
110+
std::cout << " # of clustered b hadrons: " << bHadrons.size() << std::endl;
111+
for(reco::GenParticleRefVector::const_iterator it = bHadrons.begin(); it != bHadrons.end(); ++it)
112+
{
113+
float dist = reco::deltaR( aJet->eta(), aJet->phi(), (*it)->eta(), (*it)->phi() );
114+
float dist2 = reco::deltaR( aJet->rapidity(), aJet->phi(), (*it)->rapidity(), (*it)->phi() );
115+
std::cout << " b hadron " << (it-bHadrons.begin())
116+
<< " PdgID, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId()
117+
<< ", (" << (*it)->pt()
118+
<< "," << (*it)->eta()
119+
<< "," << (*it)->rapidity()
120+
<< "," << (*it)->phi()
121+
<< "), " << dist
122+
<< ", " << dist2 << std::endl;
123+
}
124+
125+
const reco::GenParticleRefVector & cHadrons = aInfo.getcHadrons();
126+
std::cout << " # of clustered c hadrons: " << cHadrons.size() << std::endl;
127+
for(reco::GenParticleRefVector::const_iterator it = cHadrons.begin(); it != cHadrons.end(); ++it)
128+
{
129+
float dist = reco::deltaR( aJet->eta(), aJet->phi(), (*it)->eta(), (*it)->phi() );
130+
float dist2 = reco::deltaR( aJet->rapidity(), aJet->phi(), (*it)->rapidity(), (*it)->phi() );
131+
std::cout << " c hadron " << (it-cHadrons.begin())
132+
<< " PdgID, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId()
133+
<< ", (" << (*it)->pt()
134+
<< "," << (*it)->eta()
135+
<< "," << (*it)->rapidity()
136+
<< "," << (*it)->phi()
137+
<< "), " << dist
138+
<< ", " << dist2 << std::endl;
139+
}
140+
// ----------------------- Partons -------------------------------
141+
std::cout << " Parton-based flavour: " << aInfo.getPartonFlavour() << std::endl;
142+
143+
const reco::GenParticleRefVector & partons = aInfo.getPartons();
144+
std::cout << " # of clustered partons: " << partons.size() << std::endl;
145+
for(reco::GenParticleRefVector::const_iterator it = partons.begin(); it != partons.end(); ++it)
146+
{
147+
float dist = reco::deltaR( aJet->eta(), aJet->phi(), (*it)->eta(), (*it)->phi() );
148+
float dist2 = reco::deltaR( aJet->rapidity(), aJet->phi(), (*it)->rapidity(), (*it)->phi() );
149+
std::cout << " Parton " << (it-partons.begin())
150+
<< " PdgID, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId()
151+
<< ", (" << (*it)->pt()
152+
<< "," << (*it)->eta()
153+
<< "," << (*it)->rapidity()
154+
<< "," << (*it)->phi()
155+
<< "), " << dist
156+
<< ", " << dist2 << std::endl;
157+
}
158+
159+
if( useSubjets_ )
160+
{
161+
if( matchedIndices.at(j - theJetFlavourInfos->begin())<0 )
162+
{
163+
std::cout << " ----------------------- Subjet Flavour Info -----------------------" << std::endl;
164+
std::cout << " No subjets assigned to this jet" << std::endl;
165+
continue;
166+
}
167+
168+
// loop over subjets
169+
std::cout << " ----------------------- Subjet Flavour Info -----------------------" << std::endl;
170+
for(size_t s=0; s<groomedJets->at(matchedIndices.at(j - theJetFlavourInfos->begin())).numberOfDaughters(); ++s)
171+
{
172+
const edm::Ptr<reco::Candidate> & subjet = groomedJets->at(matchedIndices.at(j - theJetFlavourInfos->begin())).daughterPtr(s);
173+
174+
for ( reco::JetFlavourInfoMatchingCollection::const_iterator sj = theSubjetFlavourInfos->begin();
175+
sj != theSubjetFlavourInfos->end();
176+
++sj ) {
177+
if( subjet != edm::Ptr<reco::Candidate>((*sj).first.id(), (*sj).first.get(), (*sj).first.key()) ) continue;
178+
179+
const reco::Jet *aSubjet = (*sj).first.get();
180+
aInfo = (*sj).second;
181+
std::cout << std::setprecision(2) << std::setw(6) << std::fixed
182+
<< " [printSubjetFlavourInfo] Subjet " << s << " pt, eta, rapidity, phi, dR(eta-phi), dR(rap-phi) = "
183+
<< aSubjet->pt() << ", "
184+
<< aSubjet->eta() << ", "
185+
<< aSubjet->rapidity() << ", "
186+
<< aSubjet->phi() << ", "
187+
<< reco::deltaR( aSubjet->eta(), aSubjet->phi(), aJet->eta(), aJet->phi() ) << ", "
188+
<< reco::deltaR( aSubjet->rapidity(), aSubjet->phi(), aJet->rapidity(), aJet->phi() ) << ", "
189+
<< std::endl;
190+
// ----------------------- Hadrons -------------------------------
191+
std::cout << " Hadron-based flavour: " << aInfo.getHadronFlavour() << std::endl;
192+
193+
const reco::GenParticleRefVector & bHadrons = aInfo.getbHadrons();
194+
std::cout << " # of assigned b hadrons: " << bHadrons.size() << std::endl;
195+
for(reco::GenParticleRefVector::const_iterator it = bHadrons.begin(); it != bHadrons.end(); ++it)
196+
{
197+
float dist = reco::deltaR( aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi() );
198+
float dist2 = reco::deltaR( aSubjet->rapidity(), aSubjet->phi(), (*it)->rapidity(), (*it)->phi() );
199+
std::cout << " b hadron " << (it-bHadrons.begin())
200+
<< " PdgID, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId()
201+
<< ", (" << (*it)->pt()
202+
<< "," << (*it)->eta()
203+
<< "," << (*it)->rapidity()
204+
<< "," << (*it)->phi()
205+
<< "), " << dist
206+
<< ", " << dist2 << std::endl;
207+
}
208+
209+
const reco::GenParticleRefVector & cHadrons = aInfo.getcHadrons();
210+
std::cout << " # of assigned c hadrons: " << cHadrons.size() << std::endl;
211+
for(reco::GenParticleRefVector::const_iterator it = cHadrons.begin(); it != cHadrons.end(); ++it)
212+
{
213+
float dist = reco::deltaR( aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi() );
214+
float dist2 = reco::deltaR( aSubjet->rapidity(), aSubjet->phi(), (*it)->rapidity(), (*it)->phi() );
215+
std::cout << " c hadron " << (it-cHadrons.begin())
216+
<< " PdgID, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId()
217+
<< ", (" << (*it)->pt()
218+
<< "," << (*it)->eta()
219+
<< "," << (*it)->rapidity()
220+
<< "," << (*it)->phi()
221+
<< "), " << dist
222+
<< ", " << dist2 << std::endl;
223+
}
224+
// ----------------------- Partons -------------------------------
225+
std::cout << " Parton-based flavour: " << aInfo.getPartonFlavour() << std::endl;
226+
227+
const reco::GenParticleRefVector & partons = aInfo.getPartons();
228+
std::cout << " # of assigned partons: " << partons.size() << std::endl;
229+
for(reco::GenParticleRefVector::const_iterator it = partons.begin(); it != partons.end(); ++it)
230+
{
231+
float dist = reco::deltaR( aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi() );
232+
float dist2 = reco::deltaR( aSubjet->rapidity(), aSubjet->phi(), (*it)->rapidity(), (*it)->phi() );
233+
std::cout << " Parton " << (it-partons.begin())
234+
<< " PdgID, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId()
235+
<< ", (" << (*it)->pt()
236+
<< "," << (*it)->eta()
237+
<< "," << (*it)->rapidity()
238+
<< "," << (*it)->phi()
239+
<< "), " << dist
240+
<< ", " << dist2 << std::endl;
241+
}
242+
}
243+
}
244+
}
245+
}
246+
}
247+
248+
DEFINE_FWK_MODULE( printJetFlavourInfo );

0 commit comments

Comments
 (0)