CMS 3D CMS Logo

PseudoTopProducer.cc
Go to the documentation of this file.
2 
5 
7 #include "fastjet/ClusterSequence.hh"
8 
9 using namespace std;
10 using namespace edm;
11 using namespace reco;
12 
14  finalStateToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("finalStates"))),
15  genParticleToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("genParticles"))),
16  minLeptonPt_(pset.getParameter<double>("minLeptonPt")),
17  maxLeptonEta_(pset.getParameter<double>("maxLeptonEta")),
18  minJetPt_(pset.getParameter<double>("minJetPt")),
19  maxJetEta_(pset.getParameter<double>("maxJetEta")),
20  wMass_(pset.getParameter<double>("wMass")),
21  tMass_(pset.getParameter<double>("tMass")),
22  minLeptonPtDilepton_(pset.getParameter<double>("minLeptonPtDilepton")),
23  maxLeptonEtaDilepton_(pset.getParameter<double>("maxLeptonEtaDilepton")),
24  minDileptonMassDilepton_(pset.getParameter<double>("minDileptonMassDilepton")),
25  minLeptonPtSemilepton_(pset.getParameter<double>("minLeptonPtSemilepton")),
26  maxLeptonEtaSemilepton_(pset.getParameter<double>("maxLeptonEtaSemilepton")),
27  minVetoLeptonPtSemilepton_(pset.getParameter<double>("minVetoLeptonPtSemilepton")),
28  maxVetoLeptonEtaSemilepton_(pset.getParameter<double>("maxVetoLeptonEtaSemilepton")),
29  minMETSemiLepton_(pset.getParameter<double>("minMETSemiLepton")),
30  minMtWSemiLepton_(pset.getParameter<double>("minMtWSemiLepton"))
31 {
32  const double leptonConeSize = pset.getParameter<double>("leptonConeSize");
33  const double jetConeSize = pset.getParameter<double>("jetConeSize");
34  fjLepDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, leptonConeSize);
35  fjJetDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, jetConeSize);
36 
38 
39  produces<reco::GenParticleCollection>("neutrinos");
40  produces<reco::GenJetCollection>("leptons");
41  produces<reco::GenJetCollection>("jets");
42 
43  produces<reco::GenParticleCollection>();
44 
45 }
46 
48 {
49  edm::Handle<edm::View<reco::Candidate> > finalStateHandle;
50  event.getByToken(finalStateToken_, finalStateHandle);
51 
52  edm::Handle<edm::View<reco::Candidate> > genParticleHandle;
53  event.getByToken(genParticleToken_, genParticleHandle);
54 
55  std::unique_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
56  std::unique_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
57  std::unique_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
58  auto neutrinosRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("neutrinos");
59  auto leptonsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("leptons");
60  auto jetsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("jets");
61 
62  std::unique_ptr<reco::GenParticleCollection> pseudoTop(new reco::GenParticleCollection);
63  auto pseudoTopRefHandle = event.getRefBeforePut<reco::GenParticleCollection>();
64 
65  // Collect unstable B-hadrons
66  std::set<int> bHadronIdxs;
67  for ( size_t i=0, n=genParticleHandle->size(); i<n; ++i ) {
68  const reco::Candidate& p = genParticleHandle->at(i);
69  const int status = p.status();
70  if ( status == 1 ) continue;
71 
72  // Collect B-hadrons, to be used in b tagging
73  if ( isBHadron(&p) ) bHadronIdxs.insert(-i-1);
74  }
75 
76  // Collect stable leptons and neutrinos
77  size_t nStables = 0;
78  std::vector<size_t> leptonIdxs;
79  for ( size_t i=0, n=finalStateHandle->size(); i<n; ++i ) {
80  const reco::Candidate& p = finalStateHandle->at(i);
81  const int absPdgId = abs(p.pdgId());
82  if ( p.status() != 1 ) continue;
83 
84  ++nStables;
85  if ( p.numberOfMothers() == 0 ) continue; // Skip orphans (if exists)
86  if ( p.mother()->status() == 4 ) continue; // Treat particle as hadronic if directly from the incident beam (protect orphans in MINIAOD)
87  if ( isFromHadron(&p) ) continue;
88  switch ( absPdgId ) {
89  case 11: case 13: // Leptons
90  case 22: // Photons
91  leptonIdxs.push_back(i);
92  break;
93  case 12: case 14: case 16:
94  neutrinos->push_back(reco::GenParticle(p.charge(), p.p4(), p.vertex(), p.pdgId(), p.status(), true));
95  break;
96  }
97  }
98 
99  // Sort neutrinos by pT.
100  std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
101 
102  // Make dressed leptons with anti-kt(0.1) algorithm
104  std::vector<fastjet::PseudoJet> fjLepInputs;
105  fjLepInputs.reserve(leptonIdxs.size());
106  for ( auto index : leptonIdxs ) {
107  const reco::Candidate& p = finalStateHandle->at(index);
108  if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
109 
110  fjLepInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
111  fjLepInputs.back().set_user_index(index);
112  }
113 
115  fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *fjLepDef_);
116  std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(minLeptonPt_));
117 
119  leptons->reserve(fjLepJets.size());
120  std::set<size_t> lepDauIdxs; // keep lepton constituents to remove from GenJet construction
121  for ( auto& fjJet : fjLepJets ) {
122  if ( abs(fjJet.eta()) > maxLeptonEta_ ) continue;
123 
124  // Get jet constituents from fastJet
125  const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
126  // Convert to CandidatePtr
127  std::vector<reco::CandidatePtr> constituents;
128  reco::CandidatePtr lepCand;
129  for ( auto& fjConstituent : fjConstituents ) {
130  const size_t index = fjConstituent.user_index();
131  reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
132  const int absPdgId = abs(cand->pdgId());
133  if ( absPdgId == 11 or absPdgId == 13 ) {
134  if ( lepCand.isNonnull() and lepCand->pt() > cand->pt() ) continue; // Choose one with highest pt
135  lepCand = cand;
136  }
137  constituents.push_back(cand);
138  }
139  if ( lepCand.isNull() ) continue;
140 
141  const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
142  reco::GenJet lepJet;
143  reco::writeSpecific(lepJet, jetP4, genVertex_, constituents, eventSetup);
144 
145  lepJet.setPdgId(lepCand->pdgId());
146  lepJet.setCharge(lepCand->charge());
147 
148  const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
149  lepJet.setJetArea(jetArea);
150 
151  leptons->push_back(lepJet);
152 
153  // Keep constituent indices to be used in the next step.
154  for ( auto& fjConstituent : fjConstituents ) {
155  lepDauIdxs.insert(fjConstituent.user_index());
156  }
157  }
158 
159  // Now proceed to jets.
160  // Jets: anti-kt excluding the e, mu, nu, and photons in selected leptons.
162  std::vector<fastjet::PseudoJet> fjJetInputs;
163  fjJetInputs.reserve(nStables);
164  for ( size_t i=0, n=finalStateHandle->size(); i<n; ++i ) {
165  const reco::Candidate& p = finalStateHandle->at(i);
166  if ( p.status() != 1 ) continue;
167  if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
168 
169  const int absId = std::abs(p.pdgId());
170  if ( absId == 12 or absId == 14 or absId == 16 ) continue;
171  if ( lepDauIdxs.find(i) != lepDauIdxs.end() ) continue;
172 
173  fjJetInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
174  fjJetInputs.back().set_user_index(i);
175  }
177  for ( auto index : bHadronIdxs ) {
178  const reco::Candidate& p = genParticleHandle->at(abs(index+1));
179  if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
180 
181  const double scale = 1e-20/p.p();
182  fjJetInputs.push_back(fastjet::PseudoJet(p.px()*scale, p.py()*scale, p.pz()*scale, p.energy()*scale));
183  fjJetInputs.back().set_user_index(index);
184  }
185 
187  fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *fjJetDef_);
188  std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(minJetPt_));
189 
191  jets->reserve(fjJets.size());
192  std::vector<size_t> bjetIdxs, ljetIdxs;
193  for ( auto& fjJet : fjJets ) {
194  if ( abs(fjJet.eta()) > maxJetEta_ ) continue;
195 
196  // Get jet constituents from fastJet
197  const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
198  // Convert to CandidatePtr
199  std::vector<reco::CandidatePtr> constituents;
200  bool hasBHadron = false;
201  for ( size_t j=0, m=fjConstituents.size(); j<m; ++j ) {
202  const int index = fjConstituents[j].user_index();
203  if ( bHadronIdxs.find(index) != bHadronIdxs.end() ) {
204  hasBHadron = true;
205  continue;
206  }
207  reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
208  constituents.push_back(cand);
209  }
210 
211  const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
212  reco::GenJet genJet;
213  reco::writeSpecific(genJet, jetP4, genVertex_, constituents, eventSetup);
214 
215  const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
216  genJet.setJetArea(jetArea);
217  if ( hasBHadron ) {
218  genJet.setPdgId(5);
219  bjetIdxs.push_back(jets->size());
220  }
221  else {
222  ljetIdxs.push_back(jets->size());
223  }
224 
225  jets->push_back(genJet);
226  }
227 
228  // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
229  // NOTE : A C++ trick, use do-while instead of long-nested if-statements.
230  do {
231  if ( bjetIdxs.size() < 2 ) break;
232 
233  // Note : we will do dilepton or semilepton channel only
234  if ( leptons->size() == 2 and neutrinos->size() >= 2 ) {
235  // Start from dilepton channel
236  const int q1 = leptons->at(0).charge();
237  const int q2 = leptons->at(1).charge();
238  if ( q1*q2 > 0 ) break;
239 
240  const auto& lepton1 = q1 > 0 ? leptons->at(0) : leptons->at(1);
241  const auto& lepton2 = q1 > 0 ? leptons->at(1) : leptons->at(0);
242 
243  if ( lepton1.pt() < minLeptonPtDilepton_ or std::abs(lepton1.eta()) > maxLeptonEtaDilepton_ ) break;
244  if ( lepton2.pt() < minLeptonPtDilepton_ or std::abs(lepton2.eta()) > maxLeptonEtaDilepton_ ) break;
245  if ( (lepton1.p4()+lepton2.p4()).mass() < minDileptonMassDilepton_ ) break;
246 
247  double dm = 1e9;
248  int selNu1 = -1, selNu2 = -1;
249  for ( int i=0; i<2; ++i ) { // Consider leading 2 neutrinos. Neutrino vector is already sorted by pT
250  const double dm1 = std::abs((lepton1.p4()+neutrinos->at(i).p4()).mass()-wMass_);
251  for ( int j=0; j<2; ++j ) {
252  if ( i == j ) continue;
253  const double dm2 = std::abs((lepton2.p4()+neutrinos->at(j).p4()).mass()-wMass_);
254  const double newDm = dm1+dm2;
255 
256  if ( newDm < dm ) { dm = newDm; selNu1 = i; selNu2 = j; }
257  }
258  }
259  if ( dm >= 1e9 ) break;
260 
261  const auto& nu1 = neutrinos->at(selNu1);
262  const auto& nu2 = neutrinos->at(selNu2);
263  const auto w1LVec = lepton1.p4()+nu1.p4();
264  const auto w2LVec = lepton2.p4()+nu2.p4();
265 
266  // Contiue to top quarks
267  dm = 1e9; // Reset once again for top combination.
268  int selB1 = -1, selB2 = -1;
269  for ( auto i : bjetIdxs ) {
270  const double dm1 = std::abs((w1LVec+jets->at(i).p4()).mass()-tMass_);
271  for ( auto j : bjetIdxs ) {
272  if ( i == j ) continue;
273  const double dm2 = std::abs((w2LVec+jets->at(j).p4()).mass()-tMass_);
274  const double newDm = dm1+dm2;
275 
276  if ( newDm < dm ) { dm = newDm; selB1 = i; selB2 = j; }
277  }
278  }
279  if ( dm >= 1e9 ) break;
280 
281  const auto& bJet1 = jets->at(selB1);
282  const auto& bJet2 = jets->at(selB2);
283  const auto t1LVec = w1LVec + bJet1.p4();
284  const auto t2LVec = w2LVec + bJet2.p4();
285 
286  // Recalculate lepton charges (needed due to initialization of leptons wrt sign of a charge)
287  const int lepQ1 = lepton1.charge();
288  const int lepQ2 = lepton2.charge();
289 
290  // Put all of them into candidate collection
291  reco::GenParticle t1(lepQ1*2/3., t1LVec, genVertex_, lepQ1*6, 3, false);
292  reco::GenParticle w1(lepQ1, w1LVec, genVertex_, lepQ1*24, 3, true);
293  reco::GenParticle b1(-lepQ1/3., bJet1.p4(), genVertex_, lepQ1*5, 1, true);
294  reco::GenParticle l1(lepQ1, lepton1.p4(), genVertex_, lepton1.pdgId(), 1, true);
295  reco::GenParticle n1(0, nu1.p4(), genVertex_, nu1.pdgId(), 1, true);
296 
297  reco::GenParticle t2(lepQ2*2/3., t2LVec, genVertex_, lepQ2*6, 3, false);
298  reco::GenParticle w2(lepQ2, w2LVec, genVertex_, lepQ2*24, 3, true);
299  reco::GenParticle b2(-lepQ2/3., bJet2.p4(), genVertex_, lepQ2*5, 1, true);
300  reco::GenParticle l2(lepQ2, lepton2.p4(), genVertex_, lepton2.pdgId(), 1, true);
301  reco::GenParticle n2(0, nu2.p4(), genVertex_, nu2.pdgId(), 1, true);
302 
303  pseudoTop->push_back(t1);
304  pseudoTop->push_back(t2);
305 
306  pseudoTop->push_back(w1);
307  pseudoTop->push_back(b1);
308  pseudoTop->push_back(l1);
309  pseudoTop->push_back(n1);
310 
311  pseudoTop->push_back(w2);
312  pseudoTop->push_back(b2);
313  pseudoTop->push_back(l2);
314  pseudoTop->push_back(n2);
315  }
316  else if ( !leptons->empty() and !neutrinos->empty() and ljetIdxs.size() >= 2 ) {
317  // Then continue to the semi-leptonic channel
318  const auto& lepton = leptons->at(0);
319  if ( lepton.pt() < minLeptonPtSemilepton_ or std::abs(lepton.eta()) > maxLeptonEtaSemilepton_ ) break;
320 
321  // Skip event if there are veto leptons
322  bool hasVetoLepton = false;
323  for ( auto vetoLeptonItr = next(leptons->begin()); vetoLeptonItr != leptons->end(); ++vetoLeptonItr ) {
324  if ( vetoLeptonItr->pt() > minVetoLeptonPtSemilepton_ and
325  std::abs(vetoLeptonItr->eta()) < maxVetoLeptonEtaSemilepton_ ) {
326  hasVetoLepton = true;
327  }
328  }
329  if ( hasVetoLepton ) break;
330 
331  // Calculate MET
332  double metX = 0, metY = 0;
333  for ( auto neutrino : *neutrinos ) {
334  metX += neutrino.px();
335  metY += neutrino.py();
336  }
337  const double metPt = std::hypot(metX, metY);
338  if ( metPt < minMETSemiLepton_ ) break;
339 
340  const double mtW = std::sqrt( 2*lepton.pt()*metPt*cos(lepton.phi()-atan2(metX, metY)) );
341  if ( mtW < minMtWSemiLepton_ ) break;
342 
343  // Solve pz to give wMass_^2 = (lepE+energy)^2 - (lepPx+metX)^2 - (lepPy+metY)^2 - (lepPz+pz)^2
344  // -> (wMass_^2)/2 = lepE*sqrt(metPt^2+pz^2) - lepPx*metX - lepPy*metY - lepPz*pz
345  // -> C := (wMass_^2)/2 + lepPx*metX + lepPy*metY
346  // -> lepE^2*(metPt^2+pz^2) = C^2 + 2*C*lepPz*pz + lepPz^2*pz^2
347  // -> (lepE^2-lepPz^2)*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
348  // -> lepPt^2*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
349  const double lpz = lepton.pz(), le = lepton.energy(), lpt = lepton.pt();
350  const double cf = (wMass_*wMass_)/2 + lepton.px()*metX + lepton.py()*metY;
351  const double det = cf*cf*lpz*lpz - lpt*lpt*(le*le*metPt*metPt - cf*cf);
352  const double pz = (cf*lpz + (cf < 0 ? -sqrt(det) : sqrt(det)))/lpt/lpt;
353  const reco::Candidate::LorentzVector nu1P4(metX, metY, pz, std::hypot(metPt, pz));
354  const auto w1LVec = lepton.p4()+nu1P4;
355 
356  // Continue to build leptonic pseudo top
357  double minDR = 1e9;
358  int selB1 = -1;
359  for ( auto b1Itr = bjetIdxs.begin(); b1Itr != bjetIdxs.end(); ++b1Itr ) {
360  const double dR = deltaR(jets->at(*b1Itr).p4(), w1LVec);
361  if ( dR < minDR ) { selB1 = *b1Itr; minDR = dR; }
362  }
363  if ( selB1 == -1 ) break;
364  const auto& bJet1 = jets->at(selB1);
365  const auto t1LVec = w1LVec + bJet1.p4();
366 
367  // Build hadronic pseudo W, take leading two jets (ljetIdxs are (implicitly) sorted by pT)
368  const auto& wJet1 = jets->at(ljetIdxs[0]);
369  const auto& wJet2 = jets->at(ljetIdxs[1]);
370  const auto w2LVec = wJet1.p4() + wJet2.p4();
371 
372  // Contiue to top quarks
373  double dm = 1e9;
374  int selB2 = -1;
375  for ( auto i : bjetIdxs ) {
376  if ( int(i) == selB1 ) continue;
377  const double newDm = std::abs((w2LVec+jets->at(i).p4()).mass()-tMass_);
378  if ( newDm < dm ) { dm = newDm; selB2 = i; }
379  }
380  if ( dm >= 1e9 ) break;
381 
382  const auto& bJet2 = jets->at(selB2);
383  const auto t2LVec = w2LVec + bJet2.p4();
384 
385  const int q = lepton.charge();
386  // Put all of them into candidate collection
387  reco::GenParticle t1(q*2/3., t1LVec, genVertex_, q*6, 3, false);
388  reco::GenParticle w1(q, w1LVec, genVertex_, q*24, 3, true);
389  reco::GenParticle b1(-q/3., bJet1.p4(), genVertex_, q*5, 1, true);
390  reco::GenParticle l1(q, lepton.p4(), genVertex_, lepton.pdgId(), 1, true);
391  reco::GenParticle n1(0, nu1P4, genVertex_, q*(std::abs(lepton.pdgId())+1), 1, true);
392 
393  reco::GenParticle t2(-q*2/3., t2LVec, genVertex_, -q*6, 3, false);
394  reco::GenParticle w2(-q, w2LVec, genVertex_, -q*24, 3, true);
395  reco::GenParticle b2(0, bJet2.p4(), genVertex_, -q*5, 1, true);
396  reco::GenParticle u2(0, wJet1.p4(), genVertex_, -2*q, 1, true);
397  reco::GenParticle d2(0, wJet2.p4(), genVertex_, q, 1, true);
398 
399  pseudoTop->push_back(t1);
400  pseudoTop->push_back(t2);
401 
402  pseudoTop->push_back(w1);
403  pseudoTop->push_back(b1);
404  pseudoTop->push_back(l1);
405  pseudoTop->push_back(n1);
406 
407  pseudoTop->push_back(w2);
408  pseudoTop->push_back(b2);
409  pseudoTop->push_back(u2);
410  pseudoTop->push_back(d2);
411  }
412  } while (false);
413 
414  if ( pseudoTop->size() == 10 ) // If pseudtop decay tree is completed
415  {
416  // t->W+b
417  pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 2)); // t->W
418  pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 3)); // t->b
419  pseudoTop->at(2).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->W
420  pseudoTop->at(3).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->b
421 
422  // W->lv or W->jj
423  pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 4));
424  pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 5));
425  pseudoTop->at(4).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
426  pseudoTop->at(5).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
427 
428  // tbar->W-b
429  pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 6));
430  pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 7));
431  pseudoTop->at(6).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
432  pseudoTop->at(7).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
433 
434  // W->jj
435  pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 8));
436  pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 9));
437  pseudoTop->at(8).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
438  pseudoTop->at(9).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
439  }
440 
441  event.put(std::move(neutrinos), "neutrinos");
442  event.put(std::move(leptons), "leptons");
443  event.put(std::move(jets), "jets");
444 
445  event.put(std::move(pseudoTop));
446 }
447 
449 {
450  for ( size_t i=0, n=p->numberOfDaughters(); i<n; ++i )
451  {
452  const reco::Candidate* dau = p->daughter(i);
453  if ( p->pdgId() == dau->pdgId() ) return getLast(dau);
454  }
455  return p;
456 }
457 
459 {
460  for ( size_t i=0, n=p->numberOfMothers(); i<n; ++i )
461  {
462  const reco::Candidate* mother = p->mother(i);
463  if ( !mother ) continue;
464  if ( mother->status() == 4 or mother->numberOfMothers() == 0 ) continue; // Skip incident beam
465  const int pdgId = abs(mother->pdgId());
466 
467  if ( pdgId > 100 ) return true;
468  else if ( isFromHadron(mother) ) return true;
469  }
470  return false;
471 }
472 
474 {
475  const unsigned int absPdgId = abs(p->pdgId());
476  if ( !isBHadron(absPdgId) ) return false;
477 
478  // Do not consider this particle if it has B hadron daughter
479  // For example, B* -> B0 + photon; then we drop B* and take B0 only
480  for ( int i=0, n=p->numberOfDaughters(); i<n; ++i )
481  {
482  const reco::Candidate* dau = p->daughter(i);
483  if ( isBHadron(abs(dau->pdgId())) ) return false;
484  }
485 
486  return true;
487 }
488 
489 bool PseudoTopProducer::isBHadron(const unsigned int absPdgId) const
490 {
491  if ( absPdgId <= 100 ) return false; // Fundamental particles and MC internals
492  if ( absPdgId >= 1000000000 ) return false; // Nuclei, +-10LZZZAAAI
493 
494  // General form of PDG ID is 7 digit form
495  // +- n nr nL nq1 nq2 nq3 nJ
496  //const int nJ = absPdgId % 10; // Spin
497  const int nq3 = (absPdgId / 10) % 10;
498  const int nq2 = (absPdgId / 100) % 10;
499  const int nq1 = (absPdgId / 1000) % 10;
500 
501  if ( nq3 == 0 ) return false; // Diquarks
502  if ( nq1 == 0 and nq2 == 5 ) return true; // B mesons
503  if ( nq1 == 5 ) return true; // B baryons
504 
505  return false;
506 }
507 
509  std::auto_ptr<reco::GenParticleCollection>& outColl) const
510 {
511  reco::GenParticle pOut(*dynamic_cast<const reco::GenParticle*>(p));
512  pOut.clearMothers();
513  pOut.clearDaughters();
514  pOut.resetMothers(refHandle.id());
515  pOut.resetDaughters(refHandle.id());
516 
517  outColl->push_back(pOut);
518 
519  return reco::GenParticleRef(refHandle, outColl->size()-1);
520 }
521 
const edm::EDGetTokenT< edm::View< reco::Candidate > > genParticleToken_
bool isBHadron(const reco::Candidate *p) const
PseudoTopProducer(const edm::ParameterSet &pset)
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
T getParameter(std::string const &) const
reco::Particle::Point genVertex_
virtual double pz() const =0
z coordinate of momentum vector
const double maxLeptonEtaSemilepton_
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
const double maxLeptonEta_
const double minVetoLeptonPtSemilepton_
const double maxLeptonEtaDilepton_
std::shared_ptr< JetDef > fjJetDef_
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
double px() const final
x coordinate of momentum vector
edm::Ref< GenParticleCollection > GenParticleRef
persistent reference to a GenParticle
std::vector< GenJet > GenJetCollection
collection of GenJet objects
bool isFromHadron(const reco::Candidate *p) const
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
std::shared_ptr< JetDef > fjLepDef_
const double minMtWSemiLepton_
const double minMETSemiLepton_
reco::Particle::LorentzVector LorentzVector
double q2[4]
Definition: TauolaWrapper.h:88
void clearDaughters()
clear daughter references
const double maxJetEta_
void addDaughter(const typename daughters::value_type &)
add a daughter via a reference
void resetDaughters(const edm::ProductID &id)
set daughters product ID
const double minLeptonPtSemilepton_
virtual int status() const =0
status word
virtual double energy() const =0
energy
virtual double py() const =0
y coordinate of momentum vector
bJet1
do the razor with one or two b jets (medium CSV)
auto const T2 &decltype(t1.eta()) t2
Definition: deltaR.h:16
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
bool isnan(float x)
Definition: math.h:13
virtual int pdgId() const =0
PDG identifier.
T sqrt(T t)
Definition: SSEVec.h:18
const reco::Candidate * getLast(const reco::Candidate *p)
vector< PseudoJet > jets
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
void produce(edm::Event &event, const edm::EventSetup &eventSetup) override
bool isNull() const
Checks for null.
Definition: Ptr.h:164
math::XYZPoint Point
point in the space
Definition: Particle.h:25
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Jets made from MC generator particles.
Definition: GenJet.h:24
virtual double p() const =0
magnitude of momentum vector
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:168
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
ProductID id() const
Accessor for product ID.
Definition: RefProd.h:137
reco::GenParticleRef buildGenParticle(const reco::Candidate *p, reco::GenParticleRefProd &refHandle, std::auto_ptr< reco::GenParticleCollection > &outColl) const
virtual double pt() const =0
transverse momentum
double q1[4]
Definition: TauolaWrapper.h:87
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:37
virtual int charge() const =0
electric charge
fixed size matrix
HLT enums.
const double maxVetoLeptonEtaSemilepton_
const double minLeptonPt_
virtual const Point & vertex() const =0
vertex position
const double minDileptonMassDilepton_
virtual double px() const =0
x coordinate of momentum vector
void resetMothers(const edm::ProductID &id)
set mother product ID
const edm::EDGetTokenT< edm::View< reco::Candidate > > finalStateToken_
const double minLeptonPtDilepton_
virtual size_type numberOfDaughters() const =0
number of daughters
const double minJetPt_
void clearMothers()
clear mother references
def move(src, dest)
Definition: eostools.py:510
Definition: event.py:1
void writeSpecific(reco::CaloJet &jet, reco::Particle::LorentzVector const &p4, reco::Particle::Point const &point, std::vector< reco::CandidatePtr > const &constituents, edm::EventSetup const &c)
Definition: JetSpecific.cc:41