CMS 3D CMS Logo

GenTtbarCategorizer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TopQuarkAnalysis/TopTools
4 // Class: GenTtbarCategorizer
5 //
48 //
49 // Original Author: Johannes Hauk, Nazar Bartosik
50 // Created: Sun, 14 Jun 2015 19:42:58 GMT
51 //
52 //
53 
54 // system include files
55 #include <memory>
56 #include <algorithm>
57 #include <functional>
58 
59 // user include files
68 
69 //
70 // class declaration
71 //
72 
74 public:
75  explicit GenTtbarCategorizer(const edm::ParameterSet&);
76 
77  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
78 
79 private:
80  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
81 
82  std::vector<int> nHadronsOrderedJetIndices(const std::map<int, int>& m_jetIndex) const;
83 
84  // ----------member data ---------------------------
85 
86  // Jet configuration
87  const double genJetPtMin_;
88  const double genJetAbsEtaMax_;
89 
90  // Input tags
92 
101 
106 };
107 
108 //
109 // constants, enums and typedefs
110 //
111 
112 //
113 // static data member definitions
114 //
115 
116 //
117 // constructors and destructor
118 //
120  : genJetPtMin_(iConfig.getParameter<double>("genJetPtMin")),
121  genJetAbsEtaMax_(iConfig.getParameter<double>("genJetAbsEtaMax")),
122  genJetsToken_(consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("genJets"))),
123  genBHadJetIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadJetIndex"))),
124  genBHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFlavour"))),
125  genBHadFromTopWeakDecayToken_(
126  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFromTopWeakDecay"))),
127  genBHadPlusMothersToken_(
128  consumes<std::vector<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothers"))),
129  genBHadPlusMothersIndicesToken_(
130  consumes<std::vector<std::vector<int> > >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothersIndices"))),
131  genBHadIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadIndex"))),
132  genBHadLeptonHadronIndexToken_(
133  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadLeptonHadronIndex"))),
134  genBHadLeptonViaTauToken_(
135  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadLeptonViaTau"))),
136  genCHadJetIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadJetIndex"))),
137  genCHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadFlavour"))),
138  genCHadFromTopWeakDecayToken_(
139  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadFromTopWeakDecay"))),
140  genCHadBHadronIdToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadBHadronId"))) {
141  produces<int>("genTtbarId");
142 }
143 
144 //
145 // member functions
146 //
147 
148 // ------------ method called to produce the data ------------
150  // Access gen jets
152  iEvent.getByToken(genJetsToken_, genJets);
153 
154  // Access B hadrons information
157 
158  edm::Handle<std::vector<int> > genBHadJetIndex;
159  iEvent.getByToken(genBHadJetIndexToken_, genBHadJetIndex);
160 
163 
166 
169 
172 
173  edm::Handle<std::vector<int> > genBHadLeptonHadronIndex;
174  iEvent.getByToken(genBHadLeptonHadronIndexToken_, genBHadLeptonHadronIndex);
175 
176  edm::Handle<std::vector<int> > genBHadLeptonViaTau;
177  iEvent.getByToken(genBHadLeptonViaTauToken_, genBHadLeptonViaTau);
178 
179  // Access C hadrons information
180  edm::Handle<std::vector<int> > genCHadFlavour;
181  iEvent.getByToken(genCHadFlavourToken_, genCHadFlavour);
182 
183  edm::Handle<std::vector<int> > genCHadJetIndex;
184  iEvent.getByToken(genCHadJetIndexToken_, genCHadJetIndex);
185 
186  edm::Handle<std::vector<int> > genCHadFromTopWeakDecay;
187  iEvent.getByToken(genCHadFromTopWeakDecayToken_, genCHadFromTopWeakDecay);
188 
189  edm::Handle<std::vector<int> > genCHadBHadronId;
190  iEvent.getByToken(genCHadBHadronIdToken_, genCHadBHadronId);
191 
192  // Map <jet index, number of specific hadrons in jet>
193  // B jets with b hadrons directly from t->b decay
194  std::map<int, int> bJetFromTopIds;
195  // B jets with b hadrons from W->b decay
196  std::map<int, int> bJetFromWIds;
197  // C jets with c hadrons from W->c decay
198  std::map<int, int> cJetFromWIds;
199  // B jets with b hadrons before top quark decay chain
200  std::map<int, int> bJetAdditionalIds;
201  // C jets with c hadrons before top quark decay chain
202  std::map<int, int> cJetAdditionalIds;
203 
204  // Count number of specific b hadrons in each jet
205  for (size_t hadronId = 0; hadronId < genBHadIndex->size(); ++hadronId) {
206  // Index of jet associated to the hadron
207  const int jetIndex = genBHadJetIndex->at(hadronId);
208  // Skip hadrons which have no associated jet
209  if (jetIndex < 0)
210  continue;
211  // Skip if jet is not in acceptance
212  if (genJets->at(jetIndex).pt() < genJetPtMin_)
213  continue;
214  if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
215  continue;
216  // Flavour of the hadron's origin
217  const int flavour = genBHadFlavour->at(hadronId);
218  // Jet from t->b decay [pdgId(top)=6]
219  if (std::abs(flavour) == 6) {
220  if (bJetFromTopIds.count(jetIndex) < 1)
221  bJetFromTopIds[jetIndex] = 1;
222  else
223  bJetFromTopIds[jetIndex]++;
224  continue;
225  }
226  // Jet from W->b decay [pdgId(W)=24]
227  if (std::abs(flavour) == 24) {
228  if (bJetFromWIds.count(jetIndex) < 1)
229  bJetFromWIds[jetIndex] = 1;
230  else
231  bJetFromWIds[jetIndex]++;
232  continue;
233  }
234  // Identify jets with b hadrons not from top-quark or W-boson decay
235  if (bJetAdditionalIds.count(jetIndex) < 1)
236  bJetAdditionalIds[jetIndex] = 1;
237  else
238  bJetAdditionalIds[jetIndex]++;
239  }
240 
241  // Cleaning up b jets from W->b decays
242  for (std::map<int, int>::iterator it = bJetFromWIds.begin(); it != bJetFromWIds.end();) {
243  // Cannot be a b jet from t->b decay
244  if (bJetFromTopIds.count(it->first) > 0)
245  bJetFromWIds.erase(it++);
246  else
247  ++it;
248  }
249 
250  // Cleaning up additional b jets
251  for (std::map<int, int>::iterator it = bJetAdditionalIds.begin(); it != bJetAdditionalIds.end();) {
252  // Cannot be a b jet from t->b decay
253  if (bJetFromTopIds.count(it->first) > 0)
254  bJetAdditionalIds.erase(it++);
255  // Cannot be a b jet from W->b decay
256  else if (bJetFromWIds.count(it->first) > 0)
257  bJetAdditionalIds.erase(it++);
258  else
259  ++it;
260  }
261 
262  // Count number of specific c hadrons in each c jet
263  for (size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) {
264  // Index of jet associated to the hadron
265  const int jetIndex = genCHadJetIndex->at(hadronId);
266  // Skip hadrons which have no associated jet
267  if (jetIndex < 0)
268  continue;
269  // Skip c hadrons that are coming from b hadrons
270  if (genCHadBHadronId->at(hadronId) >= 0)
271  continue;
272  // Skip if jet is not in acceptance
273  if (genJets->at(jetIndex).pt() < genJetPtMin_)
274  continue;
275  if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
276  continue;
277  // Skip if jet contains a b hadron
278  if (bJetFromTopIds.count(jetIndex) > 0)
279  continue;
280  if (bJetFromWIds.count(jetIndex) > 0)
281  continue;
282  if (bJetAdditionalIds.count(jetIndex) > 0)
283  continue;
284  // Flavour of the hadron's origin
285  const int flavour = genCHadFlavour->at(hadronId);
286  // Jet from W->c decay [pdgId(W)=24]
287  if (std::abs(flavour) == 24) {
288  if (cJetFromWIds.count(jetIndex) < 1)
289  cJetFromWIds[jetIndex] = 1;
290  else
291  cJetFromWIds[jetIndex]++;
292  continue;
293  }
294  // Identify jets with c hadrons not from W-boson decay
295  if (cJetAdditionalIds.count(jetIndex) < 1)
296  cJetAdditionalIds[jetIndex] = 1;
297  else
298  cJetAdditionalIds[jetIndex]++;
299  }
300 
301  // Cleaning up additional c jets
302  for (std::map<int, int>::iterator it = cJetAdditionalIds.begin(); it != cJetAdditionalIds.end();) {
303  // Cannot be a c jet from W->c decay
304  if (cJetFromWIds.count(it->first) > 0)
305  cJetAdditionalIds.erase(it++);
306  else
307  ++it;
308  }
309 
310  // Categorize event based on number of additional b/c jets
311  // and number of corresponding hadrons in each of them
312  int additionalJetEventId = bJetFromTopIds.size() * 100 + bJetFromWIds.size() * 1000 + cJetFromWIds.size() * 10000;
313  // tt + 1 additional b jet
314  if (bJetAdditionalIds.size() == 1) {
315  const int nHadronsInJet = bJetAdditionalIds.begin()->second;
316  // tt + 1 additional b jet from 1 additional b hadron
317  if (nHadronsInJet == 1)
318  additionalJetEventId += 51;
319  // tt + 1 additional b jet from >=2 additional b hadrons
320  else
321  additionalJetEventId += 52;
322  }
323  // tt + >=2 additional b jets
324  else if (bJetAdditionalIds.size() > 1) {
325  // Check first two additional b jets (rare cases could have more)
326  const std::vector<int> orderedJetIndices = nHadronsOrderedJetIndices(bJetAdditionalIds);
327  int nHadronsInJet1 = bJetAdditionalIds.at(orderedJetIndices.at(0));
328  int nHadronsInJet2 = bJetAdditionalIds.at(orderedJetIndices.at(1));
329  // tt + >=2 additional b jets each from 1 additional b hadron
330  if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
331  additionalJetEventId += 53;
332  // tt + >=2 additional b jets one of which from >=2 additional b hadrons
333  else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
334  additionalJetEventId += 54;
335  // tt + >=2 additional b jets each from >=2 additional b hadrons
336  else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
337  additionalJetEventId += 55;
338  }
339  // tt + no additional b jets
340  else {
341  // tt + 1 additional c jet
342  if (cJetAdditionalIds.size() == 1) {
343  const int nHadronsInJet = cJetAdditionalIds.begin()->second;
344  // tt + 1 additional c jet from 1 additional c hadron
345  if (nHadronsInJet == 1)
346  additionalJetEventId += 41;
347  // tt + 1 additional c jet from >=2 additional c hadrons
348  else
349  additionalJetEventId += 42;
350  }
351  // tt + >=2 additional c jets
352  else if (cJetAdditionalIds.size() > 1) {
353  // Check first two additional c jets (rare cases could have more)
354  const std::vector<int> orderedJetIndices = nHadronsOrderedJetIndices(cJetAdditionalIds);
355  int nHadronsInJet1 = cJetAdditionalIds.at(orderedJetIndices.at(0));
356  int nHadronsInJet2 = cJetAdditionalIds.at(orderedJetIndices.at(1));
357  // tt + >=2 additional c jets each from 1 additional c hadron
358  if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
359  additionalJetEventId += 43;
360  // tt + >=2 additional c jets one of which from >=2 additional c hadrons
361  else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
362  additionalJetEventId += 44;
363  // tt + >=2 additional c jets each from >=2 additional c hadrons
364  else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
365  additionalJetEventId += 45;
366  }
367  // tt + no additional c jets
368  else {
369  // tt + light jets
370  additionalJetEventId += 0;
371  }
372  }
373 
374  std::unique_ptr<int> ttbarId(new int);
375  *ttbarId = additionalJetEventId;
376  iEvent.put(std::move(ttbarId), "genTtbarId");
377 }
378 
379 // ------------ method returns a vector of jet indices from the given map, sorted by N hadrons in descending order ------------
380 std::vector<int> GenTtbarCategorizer::nHadronsOrderedJetIndices(const std::map<int, int>& m_jetIndex) const {
381  const int nElements = m_jetIndex.size();
382  std::vector<std::pair<int, int> > v_jetNhadIndexPair;
383  v_jetNhadIndexPair.reserve(nElements);
384  for (std::map<int, int>::const_iterator it = m_jetIndex.begin(); it != m_jetIndex.end(); ++it) {
385  const int jetIndex = it->first;
386  const int nHadrons = it->second;
387  v_jetNhadIndexPair.push_back(std::pair<int, int>(nHadrons, jetIndex));
388  }
389  // Sorting the vector of pairs by their key value
390  std::sort(v_jetNhadIndexPair.begin(), v_jetNhadIndexPair.end(), std::greater<std::pair<int, int> >());
391  // Building the vector of indices in the proper order
392  std::vector<int> v_orderedJetIndices;
393  v_orderedJetIndices.reserve(nElements);
394  for (std::vector<std::pair<int, int> >::const_iterator it = v_jetNhadIndexPair.begin();
395  it != v_jetNhadIndexPair.end();
396  ++it) {
397  v_orderedJetIndices.push_back(it->second);
398  }
399 
400  return v_orderedJetIndices;
401 }
402 
403 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
405  //The following says we do not know what parameters are allowed so do no validation
406  // Please change this to state exactly what you do use, even if it is no parameters
408 
409  desc.add<double>("genJetPtMin", 20.);
410  desc.add<double>("genJetAbsEtaMax", 2.4);
411  desc.add<edm::InputTag>("genJets", edm::InputTag("ak4GenJets"));
412  desc.add<edm::InputTag>("genBHadJetIndex", edm::InputTag("matchGenBHadron", "genBHadJetIndex"));
413  desc.add<edm::InputTag>("genBHadFlavour", edm::InputTag("matchGenBHadron", "genBHadFlavour"));
414  desc.add<edm::InputTag>("genBHadFromTopWeakDecay", edm::InputTag("matchGenBHadron", "genBHadFromTopWeakDecay"));
415  desc.add<edm::InputTag>("genBHadPlusMothers", edm::InputTag("matchGenBHadron", "genBHadPlusMothers"));
416  desc.add<edm::InputTag>("genBHadPlusMothersIndices", edm::InputTag("matchGenBHadron", "genBHadPlusMothersIndices"));
417  desc.add<edm::InputTag>("genBHadIndex", edm::InputTag("matchGenBHadron", "genBHadIndex"));
418  desc.add<edm::InputTag>("genBHadLeptonHadronIndex", edm::InputTag("matchGenBHadron", "genBHadLeptonHadronIndex"));
419  desc.add<edm::InputTag>("genBHadLeptonViaTau", edm::InputTag("matchGenBHadron", "genBHadLeptonViaTau"));
420 
421  desc.add<edm::InputTag>("genCHadJetIndex", edm::InputTag("matchGenCHadron", "genCHadJetIndex"));
422  desc.add<edm::InputTag>("genCHadFlavour", edm::InputTag("matchGenCHadron", "genCHadFlavour"));
423  desc.add<edm::InputTag>("genCHadFromTopWeakDecay", edm::InputTag("matchGenCHadron", "genCHadFromTopWeakDecay"));
424  desc.add<edm::InputTag>("genCHadBHadronId", edm::InputTag("matchGenCHadron", "genCHadBHadronId"));
425 
426  descriptions.add("categorizeGenTtbar", desc);
427 }
428 
429 //define this as a plug-in
edm::StreamID
Definition: StreamID.h:30
GenJetCollection.h
GenHFHadronMatcher_cff.flavour
flavour
Definition: GenHFHadronMatcher_cff.py:8
GenTtbarCategorizer::genBHadPlusMothersToken_
const edm::EDGetTokenT< std::vector< reco::GenParticle > > genBHadPlusMothersToken_
Definition: GenTtbarCategorizer.cc:96
GenTtbarCategorizer::GenTtbarCategorizer
GenTtbarCategorizer(const edm::ParameterSet &)
Definition: GenTtbarCategorizer.cc:119
GenTtbarCategorizer::genBHadIndexToken_
const edm::EDGetTokenT< std::vector< int > > genBHadIndexToken_
Definition: GenTtbarCategorizer.cc:98
GenTtbarCategorizer::genCHadFlavourToken_
const edm::EDGetTokenT< std::vector< int > > genCHadFlavourToken_
Definition: GenTtbarCategorizer.cc:103
reco::GenJetCollection
std::vector< GenJet > GenJetCollection
collection of GenJet objects
Definition: GenJetCollection.h:14
min
T min(T a, T b)
Definition: MathUtil.h:58
edm::EDGetTokenT< reco::GenJetCollection >
edm
HLT enums.
Definition: AlignableModifier.h:19
GenTtbarCategorizer::genJetPtMin_
const double genJetPtMin_
Definition: GenTtbarCategorizer.cc:87
ttHFGenFilter_cfi.genBHadFlavour
genBHadFlavour
Definition: ttHFGenFilter_cfi.py:6
HLT_FULL_cff.InputTag
InputTag
Definition: HLT_FULL_cff.py:89281
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
GenTtbarCategorizer::genCHadJetIndexToken_
const edm::EDGetTokenT< std::vector< int > > genCHadJetIndexToken_
Definition: GenTtbarCategorizer.cc:102
ttHFGenFilter_cfi.genBHadFromTopWeakDecay
genBHadFromTopWeakDecay
Definition: ttHFGenFilter_cfi.py:7
GenTtbarCategorizer::genCHadBHadronIdToken_
const edm::EDGetTokenT< std::vector< int > > genCHadBHadronIdToken_
Definition: GenTtbarCategorizer.cc:105
ttHFGenFilter_cfi.genBHadPlusMothersIndices
genBHadPlusMothersIndices
Definition: ttHFGenFilter_cfi.py:9
GenTtbarCategorizer::produce
void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override
Definition: GenTtbarCategorizer.cc:149
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
GenTtbarCategorizer
Definition: GenTtbarCategorizer.cc:73
edm::Handle< reco::GenJetCollection >
GenParticle
Definition: GenParticle.py:1
GenParticle.h
MakerMacros.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
GenTtbarCategorizer::genBHadFromTopWeakDecayToken_
const edm::EDGetTokenT< std::vector< int > > genBHadFromTopWeakDecayToken_
Definition: GenTtbarCategorizer.cc:95
GenTtbarCategorizer::genBHadLeptonHadronIndexToken_
const edm::EDGetTokenT< std::vector< int > > genBHadLeptonHadronIndexToken_
Definition: GenTtbarCategorizer.cc:99
edm::global::EDProducer
Definition: EDProducer.h:32
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
edm::ParameterSet
Definition: ParameterSet.h:47
GenTtbarCategorizer::genJetAbsEtaMax_
const double genJetAbsEtaMax_
Definition: GenTtbarCategorizer.cc:88
Event.h
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
GenTtbarCategorizer::nHadronsOrderedJetIndices
std::vector< int > nHadronsOrderedJetIndices(const std::map< int, int > &m_jetIndex) const
Definition: GenTtbarCategorizer.cc:380
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
GenTtbarCategorizer::genJetsToken_
const edm::EDGetTokenT< reco::GenJetCollection > genJetsToken_
Definition: GenTtbarCategorizer.cc:91
createfilelist.int
int
Definition: createfilelist.py:10
iEvent
int iEvent
Definition: GenABIO.cc:224
GenTtbarCategorizer::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: GenTtbarCategorizer.cc:404
trackerHitRTTI::vector
Definition: trackerHitRTTI.h:21
ttbarCategorization_cff.genJets
genJets
Definition: ttbarCategorization_cff.py:29
edm::EventSetup
Definition: EventSetup.h:58
InputTag.h
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
GenTtbarCategorizer::genBHadFlavourToken_
const edm::EDGetTokenT< std::vector< int > > genBHadFlavourToken_
Definition: GenTtbarCategorizer.cc:94
Frameworkfwd.h
GenTtbarCategorizer::genBHadJetIndexToken_
const edm::EDGetTokenT< std::vector< int > > genBHadJetIndexToken_
Definition: GenTtbarCategorizer.cc:93
GenTtbarCategorizer::genBHadLeptonViaTauToken_
const edm::EDGetTokenT< std::vector< int > > genBHadLeptonViaTauToken_
Definition: GenTtbarCategorizer.cc:100
GenTtbarCategorizer::genCHadFromTopWeakDecayToken_
const edm::EDGetTokenT< std::vector< int > > genCHadFromTopWeakDecayToken_
Definition: GenTtbarCategorizer.cc:104
ttHFGenFilter_cfi.genBHadPlusMothers
genBHadPlusMothers
Definition: ttHFGenFilter_cfi.py:8
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterSet.h
EDProducer.h
edm::Event
Definition: Event.h:73
edm::InputTag
Definition: InputTag.h:15
GenTtbarCategorizer::genBHadPlusMothersIndicesToken_
const edm::EDGetTokenT< std::vector< std::vector< int > > > genBHadPlusMothersIndicesToken_
Definition: GenTtbarCategorizer.cc:97
ttHFGenFilter_cfi.genBHadIndex
genBHadIndex
Definition: ttHFGenFilter_cfi.py:10