CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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  ~GenTtbarCategorizer() override;
77 
78  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
79 
80 private:
81  void beginJob() override;
82  void produce(edm::Event&, const edm::EventSetup&) override;
83  void endJob() override;
84 
85  //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override;
86  //virtual void endRun(edm::Run const&, edm::EventSetup const&) override;
87  //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
88  //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
89 
90  std::vector<int> nHadronsOrderedJetIndices(const std::map<int, int>& m_jetIndex);
91 
92  // ----------member data ---------------------------
93 
94  // Jet configuration
95  const double genJetPtMin_;
96  const double genJetAbsEtaMax_;
97 
98  // Input tags
100 
109 
114 };
115 
116 //
117 // constants, enums and typedefs
118 //
119 
120 //
121 // static data member definitions
122 //
123 
124 //
125 // constructors and destructor
126 //
128  : genJetPtMin_(iConfig.getParameter<double>("genJetPtMin")),
129  genJetAbsEtaMax_(iConfig.getParameter<double>("genJetAbsEtaMax")),
130  genJetsToken_(consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("genJets"))),
131  genBHadJetIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadJetIndex"))),
132  genBHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFlavour"))),
134  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFromTopWeakDecay"))),
136  consumes<std::vector<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothers"))),
138  consumes<std::vector<std::vector<int> > >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothersIndices"))),
139  genBHadIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadIndex"))),
141  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadLeptonHadronIndex"))),
143  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadLeptonViaTau"))),
144  genCHadJetIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadJetIndex"))),
145  genCHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadFlavour"))),
147  consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadFromTopWeakDecay"))),
148  genCHadBHadronIdToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadBHadronId"))) {
149  produces<int>("genTtbarId");
150 }
151 
153 
154 //
155 // member functions
156 //
157 
158 // ------------ method called to produce the data ------------
160  // Access gen jets
162  iEvent.getByToken(genJetsToken_, genJets);
163 
164  // Access B hadrons information
166  iEvent.getByToken(genBHadFlavourToken_, genBHadFlavour);
167 
169  iEvent.getByToken(genBHadJetIndexToken_, genBHadJetIndex);
170 
172  iEvent.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay);
173 
175  iEvent.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers);
176 
178  iEvent.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices);
179 
181  iEvent.getByToken(genBHadIndexToken_, genBHadIndex);
182 
184  iEvent.getByToken(genBHadLeptonHadronIndexToken_, genBHadLeptonHadronIndex);
185 
187  iEvent.getByToken(genBHadLeptonViaTauToken_, genBHadLeptonViaTau);
188 
189  // Access C hadrons information
191  iEvent.getByToken(genCHadFlavourToken_, genCHadFlavour);
192 
194  iEvent.getByToken(genCHadJetIndexToken_, genCHadJetIndex);
195 
197  iEvent.getByToken(genCHadFromTopWeakDecayToken_, genCHadFromTopWeakDecay);
198 
200  iEvent.getByToken(genCHadBHadronIdToken_, genCHadBHadronId);
201 
202  // Map <jet index, number of specific hadrons in jet>
203  // B jets with b hadrons directly from t->b decay
204  std::map<int, int> bJetFromTopIds;
205  // B jets with b hadrons from W->b decay
206  std::map<int, int> bJetFromWIds;
207  // C jets with c hadrons from W->c decay
208  std::map<int, int> cJetFromWIds;
209  // B jets with b hadrons before top quark decay chain
210  std::map<int, int> bJetAdditionalIds;
211  // C jets with c hadrons before top quark decay chain
212  std::map<int, int> cJetAdditionalIds;
213 
214  // Count number of specific b hadrons in each jet
215  for (size_t hadronId = 0; hadronId < genBHadIndex->size(); ++hadronId) {
216  // Index of jet associated to the hadron
217  const int jetIndex = genBHadJetIndex->at(hadronId);
218  // Skip hadrons which have no associated jet
219  if (jetIndex < 0)
220  continue;
221  // Skip if jet is not in acceptance
222  if (genJets->at(jetIndex).pt() < genJetPtMin_)
223  continue;
224  if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
225  continue;
226  // Flavour of the hadron's origin
227  const int flavour = genBHadFlavour->at(hadronId);
228  // Jet from t->b decay [pdgId(top)=6]
229  if (std::abs(flavour) == 6) {
230  if (bJetFromTopIds.count(jetIndex) < 1)
231  bJetFromTopIds[jetIndex] = 1;
232  else
233  bJetFromTopIds[jetIndex]++;
234  continue;
235  }
236  // Jet from W->b decay [pdgId(W)=24]
237  if (std::abs(flavour) == 24) {
238  if (bJetFromWIds.count(jetIndex) < 1)
239  bJetFromWIds[jetIndex] = 1;
240  else
241  bJetFromWIds[jetIndex]++;
242  continue;
243  }
244  // Identify jets with b hadrons not from top-quark or W-boson decay
245  if (bJetAdditionalIds.count(jetIndex) < 1)
246  bJetAdditionalIds[jetIndex] = 1;
247  else
248  bJetAdditionalIds[jetIndex]++;
249  }
250 
251  // Cleaning up b jets from W->b decays
252  for (std::map<int, int>::iterator it = bJetFromWIds.begin(); it != bJetFromWIds.end();) {
253  // Cannot be a b jet from t->b decay
254  if (bJetFromTopIds.count(it->first) > 0)
255  bJetFromWIds.erase(it++);
256  else
257  ++it;
258  }
259 
260  // Cleaning up additional b jets
261  for (std::map<int, int>::iterator it = bJetAdditionalIds.begin(); it != bJetAdditionalIds.end();) {
262  // Cannot be a b jet from t->b decay
263  if (bJetFromTopIds.count(it->first) > 0)
264  bJetAdditionalIds.erase(it++);
265  // Cannot be a b jet from W->b decay
266  else if (bJetFromWIds.count(it->first) > 0)
267  bJetAdditionalIds.erase(it++);
268  else
269  ++it;
270  }
271 
272  // Count number of specific c hadrons in each c jet
273  for (size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) {
274  // Index of jet associated to the hadron
275  const int jetIndex = genCHadJetIndex->at(hadronId);
276  // Skip hadrons which have no associated jet
277  if (jetIndex < 0)
278  continue;
279  // Skip c hadrons that are coming from b hadrons
280  if (genCHadBHadronId->at(hadronId) >= 0)
281  continue;
282  // Skip if jet is not in acceptance
283  if (genJets->at(jetIndex).pt() < genJetPtMin_)
284  continue;
285  if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
286  continue;
287  // Skip if jet contains a b hadron
288  if (bJetFromTopIds.count(jetIndex) > 0)
289  continue;
290  if (bJetFromWIds.count(jetIndex) > 0)
291  continue;
292  if (bJetAdditionalIds.count(jetIndex) > 0)
293  continue;
294  // Flavour of the hadron's origin
295  const int flavour = genCHadFlavour->at(hadronId);
296  // Jet from W->c decay [pdgId(W)=24]
297  if (std::abs(flavour) == 24) {
298  if (cJetFromWIds.count(jetIndex) < 1)
299  cJetFromWIds[jetIndex] = 1;
300  else
301  cJetFromWIds[jetIndex]++;
302  continue;
303  }
304  // Identify jets with c hadrons not from W-boson decay
305  if (cJetAdditionalIds.count(jetIndex) < 1)
306  cJetAdditionalIds[jetIndex] = 1;
307  else
308  cJetAdditionalIds[jetIndex]++;
309  }
310 
311  // Cleaning up additional c jets
312  for (std::map<int, int>::iterator it = cJetAdditionalIds.begin(); it != cJetAdditionalIds.end();) {
313  // Cannot be a c jet from W->c decay
314  if (cJetFromWIds.count(it->first) > 0)
315  cJetAdditionalIds.erase(it++);
316  else
317  ++it;
318  }
319 
320  // Categorize event based on number of additional b/c jets
321  // and number of corresponding hadrons in each of them
322  int additionalJetEventId = bJetFromTopIds.size() * 100 + bJetFromWIds.size() * 1000 + cJetFromWIds.size() * 10000;
323  // tt + 1 additional b jet
324  if (bJetAdditionalIds.size() == 1) {
325  const int nHadronsInJet = bJetAdditionalIds.begin()->second;
326  // tt + 1 additional b jet from 1 additional b hadron
327  if (nHadronsInJet == 1)
328  additionalJetEventId += 51;
329  // tt + 1 additional b jet from >=2 additional b hadrons
330  else
331  additionalJetEventId += 52;
332  }
333  // tt + >=2 additional b jets
334  else if (bJetAdditionalIds.size() > 1) {
335  // Check first two additional b jets (rare cases could have more)
336  const std::vector<int> orderedJetIndices = nHadronsOrderedJetIndices(bJetAdditionalIds);
337  int nHadronsInJet1 = bJetAdditionalIds.at(orderedJetIndices.at(0));
338  int nHadronsInJet2 = bJetAdditionalIds.at(orderedJetIndices.at(1));
339  // tt + >=2 additional b jets each from 1 additional b hadron
340  if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
341  additionalJetEventId += 53;
342  // tt + >=2 additional b jets one of which from >=2 additional b hadrons
343  else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
344  additionalJetEventId += 54;
345  // tt + >=2 additional b jets each from >=2 additional b hadrons
346  else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
347  additionalJetEventId += 55;
348  }
349  // tt + no additional b jets
350  else {
351  // tt + 1 additional c jet
352  if (cJetAdditionalIds.size() == 1) {
353  const int nHadronsInJet = cJetAdditionalIds.begin()->second;
354  // tt + 1 additional c jet from 1 additional c hadron
355  if (nHadronsInJet == 1)
356  additionalJetEventId += 41;
357  // tt + 1 additional c jet from >=2 additional c hadrons
358  else
359  additionalJetEventId += 42;
360  }
361  // tt + >=2 additional c jets
362  else if (cJetAdditionalIds.size() > 1) {
363  // Check first two additional c jets (rare cases could have more)
364  const std::vector<int> orderedJetIndices = nHadronsOrderedJetIndices(cJetAdditionalIds);
365  int nHadronsInJet1 = cJetAdditionalIds.at(orderedJetIndices.at(0));
366  int nHadronsInJet2 = cJetAdditionalIds.at(orderedJetIndices.at(1));
367  // tt + >=2 additional c jets each from 1 additional c hadron
368  if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
369  additionalJetEventId += 43;
370  // tt + >=2 additional c jets one of which from >=2 additional c hadrons
371  else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
372  additionalJetEventId += 44;
373  // tt + >=2 additional c jets each from >=2 additional c hadrons
374  else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
375  additionalJetEventId += 45;
376  }
377  // tt + no additional c jets
378  else {
379  // tt + light jets
380  additionalJetEventId += 0;
381  }
382  }
383 
384  std::unique_ptr<int> ttbarId(new int);
385  *ttbarId = additionalJetEventId;
386  iEvent.put(std::move(ttbarId), "genTtbarId");
387 }
388 
389 // ------------ method called once each job just before starting event loop ------------
391 
392 // ------------ method called once each job just after ending the event loop ------------
394 
395 // ------------ method called when starting to processes a run ------------
396 /*
397 void
398 GenTtbarCategorizer::beginRun(edm::Run const&, edm::EventSetup const&)
399 {
400 }
401 */
402 
403 // ------------ method called when ending the processing of a run ------------
404 /*
405 void
406 GenTtbarCategorizer::endRun(edm::Run const&, edm::EventSetup const&)
407 {
408 }
409 */
410 
411 // ------------ method called when starting to processes a luminosity block ------------
412 /*
413 void
414 GenTtbarCategorizer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
415 {
416 }
417 */
418 
419 // ------------ method called when ending the processing of a luminosity block ------------
420 /*
421 void
422 GenTtbarCategorizer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
423 {
424 }
425 */
426 
427 // ------------ method returns a vector of jet indices from the given map, sorted by N hadrons in descending order ------------
428 std::vector<int> GenTtbarCategorizer::nHadronsOrderedJetIndices(const std::map<int, int>& m_jetIndex) {
429  const int nElements = m_jetIndex.size();
430  std::vector<std::pair<int, int> > v_jetNhadIndexPair;
431  v_jetNhadIndexPair.reserve(nElements);
432  for (std::map<int, int>::const_iterator it = m_jetIndex.begin(); it != m_jetIndex.end(); ++it) {
433  const int jetIndex = it->first;
434  const int nHadrons = it->second;
435  v_jetNhadIndexPair.push_back(std::pair<int, int>(nHadrons, jetIndex));
436  }
437  // Sorting the vector of pairs by their key value
438  std::sort(v_jetNhadIndexPair.begin(), v_jetNhadIndexPair.end(), std::greater<std::pair<int, int> >());
439  // Building the vector of indices in the proper order
440  std::vector<int> v_orderedJetIndices;
441  v_orderedJetIndices.reserve(nElements);
442  for (std::vector<std::pair<int, int> >::const_iterator it = v_jetNhadIndexPair.begin();
443  it != v_jetNhadIndexPair.end();
444  ++it) {
445  v_orderedJetIndices.push_back(it->second);
446  }
447 
448  return v_orderedJetIndices;
449 }
450 
451 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
453  //The following says we do not know what parameters are allowed so do no validation
454  // Please change this to state exactly what you do use, even if it is no parameters
456  desc.setUnknown();
457  descriptions.addDefault(desc);
458 }
459 
460 //define this as a plug-in
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:131
const edm::EDGetTokenT< std::vector< int > > genBHadFromTopWeakDecayToken_
const edm::EDGetTokenT< std::vector< int > > genBHadJetIndexToken_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:525
const edm::EDGetTokenT< std::vector< int > > genCHadFromTopWeakDecayToken_
std::vector< GenJet > GenJetCollection
collection of GenJet objects
const edm::EDGetTokenT< std::vector< int > > genCHadFlavourToken_
std::vector< int > nHadronsOrderedJetIndices(const std::map< int, int > &m_jetIndex)
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
int iEvent
Definition: GenABIO.cc:224
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void addDefault(ParameterSetDescription const &psetDescription)
const edm::EDGetTokenT< std::vector< int > > genBHadFlavourToken_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
T min(T a, T b)
Definition: MathUtil.h:58
const edm::EDGetTokenT< std::vector< std::vector< int > > > genBHadPlusMothersIndicesToken_
const edm::EDGetTokenT< std::vector< int > > genBHadLeptonHadronIndexToken_
const edm::EDGetTokenT< std::vector< int > > genBHadIndexToken_
const edm::EDGetTokenT< std::vector< int > > genBHadLeptonViaTauToken_
fixed size matrix
HLT enums.
const edm::EDGetTokenT< std::vector< reco::GenParticle > > genBHadPlusMothersToken_
GenTtbarCategorizer(const edm::ParameterSet &)
void produce(edm::Event &, const edm::EventSetup &) override
const edm::EDGetTokenT< std::vector< int > > genCHadJetIndexToken_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const edm::EDGetTokenT< std::vector< int > > genCHadBHadronIdToken_
const edm::EDGetTokenT< reco::GenJetCollection > genJetsToken_
def move(src, dest)
Definition: eostools.py:511