CMS 3D CMS Logo
List of all members | Public Member Functions | Private Types | Private Member Functions | Private Attributes
PseudoTopProducer Class Reference

#include <PseudoTopProducer.h>

Inheritance diagram for PseudoTopProducer:
edm::stream::EDProducer<>

Public Member Functions

void produce (edm::Event &event, const edm::EventSetup &eventSetup) override
 
 PseudoTopProducer (const edm::ParameterSet &pset)
 
- Public Member Functions inherited from edm::stream::EDProducer<>
 EDProducer ()=default
 
bool hasAbilityToProduceInBeginLumis () const final
 
bool hasAbilityToProduceInBeginProcessBlocks () const final
 
bool hasAbilityToProduceInBeginRuns () const final
 
bool hasAbilityToProduceInEndLumis () const final
 
bool hasAbilityToProduceInEndProcessBlocks () const final
 
bool hasAbilityToProduceInEndRuns () const final
 

Private Types

typedef fastjet::JetDefinition JetDef
 
typedef reco::Particle::LorentzVector LorentzVector
 

Private Member Functions

reco::GenParticleRef buildGenParticle (const reco::Candidate *p, reco::GenParticleRefProd &refHandle, std::auto_ptr< reco::GenParticleCollection > &outColl) const
 
const reco::Candidate * getLast (const reco::Candidate *p)
 
void insertAllDaughters (const reco::Candidate *p, std::set< const reco::Candidate * > &list) const
 
bool isBHadron (const reco::Candidate *p) const
 
bool isBHadron (const unsigned int pdgId) const
 
bool isFromHadron (const reco::Candidate *p) const
 

Private Attributes

const edm::EDGetTokenT< edm::View< reco::Candidate > > finalStateToken_
 
std::shared_ptr< JetDef > fjJetDef_
 
std::shared_ptr< JetDef > fjLepDef_
 
const edm::EDGetTokenT< edm::View< reco::Candidate > > genParticleToken_
 
reco::Particle::Point genVertex_
 
const double maxJetEta_
 
const double maxLeptonEta_
 
const double maxLeptonEtaDilepton_
 
const double maxLeptonEtaSemilepton_
 
const double maxVetoLeptonEtaSemilepton_
 
const double minDileptonMassDilepton_
 
const double minJetPt_
 
const double minLeptonPt_
 
const double minLeptonPtDilepton_
 
const double minLeptonPtSemilepton_
 
const double minMETSemiLepton_
 
const double minMtWSemiLepton_
 
const double minVetoLeptonPtSemilepton_
 
const double tMass_
 
const double wMass_
 

Additional Inherited Members

- Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T... > CacheTypes
 
typedef CacheTypes::GlobalCache GlobalCache
 
typedef AbilityChecker< T... > HasAbility
 
typedef CacheTypes::LuminosityBlockCache LuminosityBlockCache
 
typedef LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache > LuminosityBlockContext
 
typedef CacheTypes::LuminosityBlockSummaryCache LuminosityBlockSummaryCache
 
typedef CacheTypes::RunCache RunCache
 
typedef RunContextT< RunCache, GlobalCache > RunContext
 
typedef CacheTypes::RunSummaryCache RunSummaryCache
 

Detailed Description

Definition at line 15 of file PseudoTopProducer.h.

Member Typedef Documentation

◆ JetDef

typedef fastjet::JetDefinition PseudoTopProducer::JetDef
private

Definition at line 43 of file PseudoTopProducer.h.

◆ LorentzVector

typedef reco::Particle::LorentzVector PseudoTopProducer::LorentzVector
private

Definition at line 30 of file PseudoTopProducer.h.

Constructor & Destructor Documentation

◆ PseudoTopProducer()

PseudoTopProducer::PseudoTopProducer ( const edm::ParameterSet &  pset)

Definition at line 13 of file PseudoTopProducer.cc.

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  const double leptonConeSize = pset.getParameter<double>("leptonConeSize");
32  const double jetConeSize = pset.getParameter<double>("jetConeSize");
33  fjLepDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, leptonConeSize);
34  fjJetDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, jetConeSize);
35 
36  genVertex_ = reco::Particle::Point(0, 0, 0);
37 
38  produces<reco::GenParticleCollection>("neutrinos");
39  produces<reco::GenJetCollection>("leptons");
40  produces<reco::GenJetCollection>("jets");
41 
42  produces<reco::GenParticleCollection>();
43 }

References fjJetDef_, fjLepDef_, genVertex_, particleLevel_cfi::jetConeSize, pseudoTop_cfi::leptonConeSize, and muonDTDigis_cfi::pset.

Member Function Documentation

◆ buildGenParticle()

reco::GenParticleRef PseudoTopProducer::buildGenParticle ( const reco::Candidate *  p,
reco::GenParticleRefProd &  refHandle,
std::auto_ptr< reco::GenParticleCollection > &  outColl 
) const
private

Definition at line 557 of file PseudoTopProducer.cc.

559  {
560  reco::GenParticle pOut(*dynamic_cast<const reco::GenParticle*>(p));
561  pOut.clearMothers();
562  pOut.clearDaughters();
563  pOut.resetMothers(refHandle.id());
564  pOut.resetDaughters(refHandle.id());
565 
566  outColl->push_back(pOut);
567 
568  return reco::GenParticleRef(refHandle, outColl->size() - 1);
569 }

References reco::CompositeRefCandidateT< D >::clearDaughters(), reco::CompositeRefCandidateT< D >::clearMothers(), edm::RefProd< C >::id(), ecalCompactTrigPrim_cfi::outColl, AlCaHLTBitMon_ParallelJobs::p, reco::CompositeRefCandidateT< D >::resetDaughters(), and reco::CompositeRefCandidateT< D >::resetMothers().

◆ getLast()

const reco::Candidate * PseudoTopProducer::getLast ( const reco::Candidate *  p)
private

Definition at line 492 of file PseudoTopProducer.cc.

492  {
493  for (size_t i = 0, n = p->numberOfDaughters(); i < n; ++i) {
494  const reco::Candidate* dau = p->daughter(i);
495  if (p->pdgId() == dau->pdgId())
496  return getLast(dau);
497  }
498  return p;
499 }

References mps_fire::i, dqmiodumpmetadata::n, AlCaHLTBitMon_ParallelJobs::p, and reco::Candidate::pdgId().

◆ insertAllDaughters()

void PseudoTopProducer::insertAllDaughters ( const reco::Candidate *  p,
std::set< const reco::Candidate * > &  list 
) const
private

◆ isBHadron() [1/2]

bool PseudoTopProducer::isBHadron ( const reco::Candidate *  p) const
private

Definition at line 518 of file PseudoTopProducer.cc.

518  {
519  const unsigned int absPdgId = abs(p->pdgId());
520  if (!isBHadron(absPdgId))
521  return false;
522 
523  // Do not consider this particle if it has B hadron daughter
524  // For example, B* -> B0 + photon; then we drop B* and take B0 only
525  for (int i = 0, n = p->numberOfDaughters(); i < n; ++i) {
526  const reco::Candidate* dau = p->daughter(i);
527  if (isBHadron(abs(dau->pdgId())))
528  return false;
529  }
530 
531  return true;
532 }

References funct::abs(), mps_fire::i, dqmiodumpmetadata::n, AlCaHLTBitMon_ParallelJobs::p, and reco::Candidate::pdgId().

Referenced by produce().

◆ isBHadron() [2/2]

bool PseudoTopProducer::isBHadron ( const unsigned int  pdgId) const
private

Definition at line 534 of file PseudoTopProducer.cc.

534  {
535  if (absPdgId <= 100)
536  return false; // Fundamental particles and MC internals
537  if (absPdgId >= 1000000000)
538  return false; // Nuclei, +-10LZZZAAAI
539 
540  // General form of PDG ID is 7 digit form
541  // +- n nr nL nq1 nq2 nq3 nJ
542  //const int nJ = absPdgId % 10; // Spin
543  const int nq3 = (absPdgId / 10) % 10;
544  const int nq2 = (absPdgId / 100) % 10;
545  const int nq1 = (absPdgId / 1000) % 10;
546 
547  if (nq3 == 0)
548  return false; // Diquarks
549  if (nq1 == 0 and nq2 == 5)
550  return true; // B mesons
551  if (nq1 == 5)
552  return true; // B baryons
553 
554  return false;
555 }

◆ isFromHadron()

bool PseudoTopProducer::isFromHadron ( const reco::Candidate *  p) const
private

Definition at line 501 of file PseudoTopProducer.cc.

501  {
502  for (size_t i = 0, n = p->numberOfMothers(); i < n; ++i) {
503  const reco::Candidate* mother = p->mother(i);
504  if (!mother)
505  continue;
506  if (mother->status() == 4 or mother->numberOfMothers() == 0)
507  continue; // Skip incident beam
508  const int pdgId = abs(mother->pdgId());
509 
510  if (pdgId > 100)
511  return true;
512  else if (isFromHadron(mother))
513  return true;
514  }
515  return false;
516 }

References funct::abs(), mps_fire::i, dqmiodumpmetadata::n, reco::Candidate::numberOfMothers(), or, AlCaHLTBitMon_ParallelJobs::p, EgammaValidation_cff::pdgId, reco::Candidate::pdgId(), and reco::Candidate::status().

Referenced by produce().

◆ produce()

void PseudoTopProducer::produce ( edm::Event &  event,
const edm::EventSetup &  eventSetup 
)
override

Build jets

Definition at line 45 of file PseudoTopProducer.cc.

45  {
46  edm::Handle<edm::View<reco::Candidate> > finalStateHandle;
47  event.getByToken(finalStateToken_, finalStateHandle);
48 
49  edm::Handle<edm::View<reco::Candidate> > genParticleHandle;
50  event.getByToken(genParticleToken_, genParticleHandle);
51 
52  std::unique_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
53  std::unique_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
54  std::unique_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
55  auto neutrinosRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("neutrinos");
56  auto leptonsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("leptons");
57  auto jetsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("jets");
58 
59  std::unique_ptr<reco::GenParticleCollection> pseudoTop(new reco::GenParticleCollection);
60  auto pseudoTopRefHandle = event.getRefBeforePut<reco::GenParticleCollection>();
61 
62  // Collect unstable B-hadrons
63  std::set<int> bHadronIdxs;
64  for (size_t i = 0, n = genParticleHandle->size(); i < n; ++i) {
65  const reco::Candidate& p = genParticleHandle->at(i);
66  const int status = p.status();
67  if (status == 1)
68  continue;
69 
70  // Collect B-hadrons, to be used in b tagging
71  if (isBHadron(&p))
72  bHadronIdxs.insert(-i - 1);
73  }
74 
75  // Collect stable leptons and neutrinos
76  size_t nStables = 0;
77  std::vector<size_t> leptonIdxs;
78  for (size_t i = 0, n = finalStateHandle->size(); i < n; ++i) {
79  const reco::Candidate& p = finalStateHandle->at(i);
80  const int absPdgId = abs(p.pdgId());
81  if (p.status() != 1)
82  continue;
83 
84  ++nStables;
85  if (p.numberOfMothers() == 0)
86  continue; // Skip orphans (if exists)
87  if (p.mother()->status() == 4)
88  continue; // Treat particle as hadronic if directly from the incident beam (protect orphans in MINIAOD)
89  if (isFromHadron(&p))
90  continue;
91  switch (absPdgId) {
92  case 11:
93  case 13: // Leptons
94  case 22: // Photons
95  leptonIdxs.push_back(i);
96  break;
97  case 12:
98  case 14:
99  case 16:
100  neutrinos->push_back(reco::GenParticle(p.charge(), p.p4(), p.vertex(), p.pdgId(), p.status(), true));
101  break;
102  }
103  }
104 
105  // Sort neutrinos by pT.
106  std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
107 
108  // Make dressed leptons with anti-kt(0.1) algorithm
110  std::vector<fastjet::PseudoJet> fjLepInputs;
111  fjLepInputs.reserve(leptonIdxs.size());
112  for (auto index : leptonIdxs) {
113  const reco::Candidate& p = finalStateHandle->at(index);
114  if (std::isnan(p.pt()) or p.pt() <= 0)
115  continue;
116 
117  fjLepInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
118  fjLepInputs.back().set_user_index(index);
119  }
120 
122  fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *fjLepDef_);
123  std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(minLeptonPt_));
124 
126  leptons->reserve(fjLepJets.size());
127  std::set<size_t> lepDauIdxs; // keep lepton constituents to remove from GenJet construction
128  for (auto& fjJet : fjLepJets) {
129  if (abs(fjJet.eta()) > maxLeptonEta_)
130  continue;
131 
132  // Get jet constituents from fastJet
133  const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
134  // Convert to CandidatePtr
135  std::vector<reco::CandidatePtr> constituents;
136  reco::CandidatePtr lepCand;
137  for (auto& fjConstituent : fjConstituents) {
138  const size_t index = fjConstituent.user_index();
139  reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
140  const int absPdgId = abs(cand->pdgId());
141  if (absPdgId == 11 or absPdgId == 13) {
142  if (lepCand.isNonnull() and lepCand->pt() > cand->pt())
143  continue; // Choose one with highest pt
144  lepCand = cand;
145  }
146  constituents.push_back(cand);
147  }
148  if (lepCand.isNull())
149  continue;
150 
151  const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
152  reco::GenJet lepJet;
153  reco::writeSpecific(lepJet, jetP4, genVertex_, constituents, eventSetup);
154 
155  lepJet.setPdgId(lepCand->pdgId());
156  lepJet.setCharge(lepCand->charge());
157 
158  const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
159  lepJet.setJetArea(jetArea);
160 
161  leptons->push_back(lepJet);
162 
163  // Keep constituent indices to be used in the next step.
164  for (auto& fjConstituent : fjConstituents) {
165  lepDauIdxs.insert(fjConstituent.user_index());
166  }
167  }
168 
169  // Now proceed to jets.
170  // Jets: anti-kt excluding the e, mu, nu, and photons in selected leptons.
172  std::vector<fastjet::PseudoJet> fjJetInputs;
173  fjJetInputs.reserve(nStables);
174  for (size_t i = 0, n = finalStateHandle->size(); i < n; ++i) {
175  const reco::Candidate& p = finalStateHandle->at(i);
176  if (p.status() != 1)
177  continue;
178  if (std::isnan(p.pt()) or p.pt() <= 0)
179  continue;
180 
181  const int absId = std::abs(p.pdgId());
182  if (absId == 12 or absId == 14 or absId == 16)
183  continue;
184  if (lepDauIdxs.find(i) != lepDauIdxs.end())
185  continue;
186 
187  fjJetInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
188  fjJetInputs.back().set_user_index(i);
189  }
191  for (auto index : bHadronIdxs) {
192  const reco::Candidate& p = genParticleHandle->at(abs(index + 1));
193  if (std::isnan(p.pt()) or p.pt() <= 0)
194  continue;
195 
196  const double scale = 1e-20 / p.p();
197  fjJetInputs.push_back(fastjet::PseudoJet(p.px() * scale, p.py() * scale, p.pz() * scale, p.energy() * scale));
198  fjJetInputs.back().set_user_index(index);
199  }
200 
202  fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *fjJetDef_);
203  std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(minJetPt_));
204 
206  jets->reserve(fjJets.size());
207  std::vector<size_t> bjetIdxs, ljetIdxs;
208  for (auto& fjJet : fjJets) {
209  if (abs(fjJet.eta()) > maxJetEta_)
210  continue;
211 
212  // Get jet constituents from fastJet
213  const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
214  // Convert to CandidatePtr
215  std::vector<reco::CandidatePtr> constituents;
216  bool hasBHadron = false;
217  for (size_t j = 0, m = fjConstituents.size(); j < m; ++j) {
218  const int index = fjConstituents[j].user_index();
219  if (bHadronIdxs.find(index) != bHadronIdxs.end()) {
220  hasBHadron = true;
221  continue;
222  }
223  reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
224  constituents.push_back(cand);
225  }
226 
227  const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
228  reco::GenJet genJet;
229  reco::writeSpecific(genJet, jetP4, genVertex_, constituents, eventSetup);
230 
231  const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
232  genJet.setJetArea(jetArea);
233  if (hasBHadron) {
234  genJet.setPdgId(5);
235  bjetIdxs.push_back(jets->size());
236  } else {
237  ljetIdxs.push_back(jets->size());
238  }
239 
240  jets->push_back(genJet);
241  }
242 
243  // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
244  // NOTE : A C++ trick, use do-while instead of long-nested if-statements.
245  do {
246  if (bjetIdxs.size() < 2)
247  break;
248 
249  // Note : we will do dilepton or semilepton channel only
250  if (leptons->size() == 2 and neutrinos->size() >= 2) {
251  // Start from dilepton channel
252  const int q1 = leptons->at(0).charge();
253  const int q2 = leptons->at(1).charge();
254  if (q1 * q2 > 0)
255  break;
256 
257  const auto& lepton1 = q1 > 0 ? leptons->at(0) : leptons->at(1);
258  const auto& lepton2 = q1 > 0 ? leptons->at(1) : leptons->at(0);
259 
260  if (lepton1.pt() < minLeptonPtDilepton_ or std::abs(lepton1.eta()) > maxLeptonEtaDilepton_)
261  break;
262  if (lepton2.pt() < minLeptonPtDilepton_ or std::abs(lepton2.eta()) > maxLeptonEtaDilepton_)
263  break;
264  if ((lepton1.p4() + lepton2.p4()).mass() < minDileptonMassDilepton_)
265  break;
266 
267  double dm = 1e9;
268  int selNu1 = -1, selNu2 = -1;
269  for (int i = 0; i < 2; ++i) { // Consider leading 2 neutrinos. Neutrino vector is already sorted by pT
270  const double dm1 = std::abs((lepton1.p4() + neutrinos->at(i).p4()).mass() - wMass_);
271  for (int j = 0; j < 2; ++j) {
272  if (i == j)
273  continue;
274  const double dm2 = std::abs((lepton2.p4() + neutrinos->at(j).p4()).mass() - wMass_);
275  const double newDm = dm1 + dm2;
276 
277  if (newDm < dm) {
278  dm = newDm;
279  selNu1 = i;
280  selNu2 = j;
281  }
282  }
283  }
284  if (dm >= 1e9)
285  break;
286 
287  const auto& nu1 = neutrinos->at(selNu1);
288  const auto& nu2 = neutrinos->at(selNu2);
289  const auto w1LVec = lepton1.p4() + nu1.p4();
290  const auto w2LVec = lepton2.p4() + nu2.p4();
291 
292  // Contiue to top quarks
293  dm = 1e9; // Reset once again for top combination.
294  int selB1 = -1, selB2 = -1;
295  for (auto i : bjetIdxs) {
296  const double dm1 = std::abs((w1LVec + jets->at(i).p4()).mass() - tMass_);
297  for (auto j : bjetIdxs) {
298  if (i == j)
299  continue;
300  const double dm2 = std::abs((w2LVec + jets->at(j).p4()).mass() - tMass_);
301  const double newDm = dm1 + dm2;
302 
303  if (newDm < dm) {
304  dm = newDm;
305  selB1 = i;
306  selB2 = j;
307  }
308  }
309  }
310  if (dm >= 1e9)
311  break;
312 
313  const auto& bJet1 = jets->at(selB1);
314  const auto& bJet2 = jets->at(selB2);
315  const auto t1LVec = w1LVec + bJet1.p4();
316  const auto t2LVec = w2LVec + bJet2.p4();
317 
318  // Recalculate lepton charges (needed due to initialization of leptons wrt sign of a charge)
319  const int lepQ1 = lepton1.charge();
320  const int lepQ2 = lepton2.charge();
321 
322  // Put all of them into candidate collection
323  reco::GenParticle t1(lepQ1 * 2 / 3., t1LVec, genVertex_, lepQ1 * 6, 3, false);
324  reco::GenParticle w1(lepQ1, w1LVec, genVertex_, lepQ1 * 24, 3, true);
325  reco::GenParticle b1(-lepQ1 / 3., bJet1.p4(), genVertex_, lepQ1 * 5, 1, true);
326  reco::GenParticle l1(lepQ1, lepton1.p4(), genVertex_, lepton1.pdgId(), 1, true);
327  reco::GenParticle n1(0, nu1.p4(), genVertex_, nu1.pdgId(), 1, true);
328 
329  reco::GenParticle t2(lepQ2 * 2 / 3., t2LVec, genVertex_, lepQ2 * 6, 3, false);
330  reco::GenParticle w2(lepQ2, w2LVec, genVertex_, lepQ2 * 24, 3, true);
331  reco::GenParticle b2(-lepQ2 / 3., bJet2.p4(), genVertex_, lepQ2 * 5, 1, true);
332  reco::GenParticle l2(lepQ2, lepton2.p4(), genVertex_, lepton2.pdgId(), 1, true);
333  reco::GenParticle n2(0, nu2.p4(), genVertex_, nu2.pdgId(), 1, true);
334 
335  pseudoTop->push_back(t1);
336  pseudoTop->push_back(t2);
337 
338  pseudoTop->push_back(w1);
339  pseudoTop->push_back(b1);
340  pseudoTop->push_back(l1);
341  pseudoTop->push_back(n1);
342 
343  pseudoTop->push_back(w2);
344  pseudoTop->push_back(b2);
345  pseudoTop->push_back(l2);
346  pseudoTop->push_back(n2);
347  } else if (!leptons->empty() and !neutrinos->empty() and ljetIdxs.size() >= 2) {
348  // Then continue to the semi-leptonic channel
349  const auto& lepton = leptons->at(0);
350  if (lepton.pt() < minLeptonPtSemilepton_ or std::abs(lepton.eta()) > maxLeptonEtaSemilepton_)
351  break;
352 
353  // Skip event if there are veto leptons
354  bool hasVetoLepton = false;
355  for (auto vetoLeptonItr = next(leptons->begin()); vetoLeptonItr != leptons->end(); ++vetoLeptonItr) {
356  if (vetoLeptonItr->pt() > minVetoLeptonPtSemilepton_ and
357  std::abs(vetoLeptonItr->eta()) < maxVetoLeptonEtaSemilepton_) {
358  hasVetoLepton = true;
359  }
360  }
361  if (hasVetoLepton)
362  break;
363 
364  // Calculate MET
365  double metX = 0, metY = 0;
366  for (const auto& neutrino : *neutrinos) {
367  metX += neutrino.px();
368  metY += neutrino.py();
369  }
370  const double metPt = std::hypot(metX, metY);
371  if (metPt < minMETSemiLepton_)
372  break;
373 
374  const double mtW = std::sqrt(2 * lepton.pt() * metPt * cos(lepton.phi() - atan2(metX, metY)));
375  if (mtW < minMtWSemiLepton_)
376  break;
377 
378  // Solve pz to give wMass_^2 = (lepE+energy)^2 - (lepPx+metX)^2 - (lepPy+metY)^2 - (lepPz+pz)^2
379  // -> (wMass_^2)/2 = lepE*sqrt(metPt^2+pz^2) - lepPx*metX - lepPy*metY - lepPz*pz
380  // -> C := (wMass_^2)/2 + lepPx*metX + lepPy*metY
381  // -> lepE^2*(metPt^2+pz^2) = C^2 + 2*C*lepPz*pz + lepPz^2*pz^2
382  // -> (lepE^2-lepPz^2)*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
383  // -> lepPt^2*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
384  const double lpz = lepton.pz(), le = lepton.energy(), lpt = lepton.pt();
385  const double cf = (wMass_ * wMass_) / 2 + lepton.px() * metX + lepton.py() * metY;
386  const double det = cf * cf * lpz * lpz - lpt * lpt * (le * le * metPt * metPt - cf * cf);
387  const double pz = (cf * lpz + (cf < 0 ? -sqrt(det) : sqrt(det))) / lpt / lpt;
388  const reco::Candidate::LorentzVector nu1P4(metX, metY, pz, std::hypot(metPt, pz));
389  const auto w1LVec = lepton.p4() + nu1P4;
390 
391  // Continue to build leptonic pseudo top
392  double minDR = 1e9;
393  int selB1 = -1;
394  for (auto b1Itr = bjetIdxs.begin(); b1Itr != bjetIdxs.end(); ++b1Itr) {
395  const double dR = deltaR(jets->at(*b1Itr).p4(), w1LVec);
396  if (dR < minDR) {
397  selB1 = *b1Itr;
398  minDR = dR;
399  }
400  }
401  if (selB1 == -1)
402  break;
403  const auto& bJet1 = jets->at(selB1);
404  const auto t1LVec = w1LVec + bJet1.p4();
405 
406  // Build hadronic pseudo W, take leading two jets (ljetIdxs are (implicitly) sorted by pT)
407  const auto& wJet1 = jets->at(ljetIdxs[0]);
408  const auto& wJet2 = jets->at(ljetIdxs[1]);
409  const auto w2LVec = wJet1.p4() + wJet2.p4();
410 
411  // Contiue to top quarks
412  double dm = 1e9;
413  int selB2 = -1;
414  for (auto i : bjetIdxs) {
415  if (int(i) == selB1)
416  continue;
417  const double newDm = std::abs((w2LVec + jets->at(i).p4()).mass() - tMass_);
418  if (newDm < dm) {
419  dm = newDm;
420  selB2 = i;
421  }
422  }
423  if (dm >= 1e9)
424  break;
425 
426  const auto& bJet2 = jets->at(selB2);
427  const auto t2LVec = w2LVec + bJet2.p4();
428 
429  const int q = lepton.charge();
430  // Put all of them into candidate collection
431  reco::GenParticle t1(q * 2 / 3., t1LVec, genVertex_, q * 6, 3, false);
432  reco::GenParticle w1(q, w1LVec, genVertex_, q * 24, 3, true);
433  reco::GenParticle b1(-q / 3., bJet1.p4(), genVertex_, q * 5, 1, true);
434  reco::GenParticle l1(q, lepton.p4(), genVertex_, lepton.pdgId(), 1, true);
435  reco::GenParticle n1(0, nu1P4, genVertex_, q * (std::abs(lepton.pdgId()) + 1), 1, true);
436 
437  reco::GenParticle t2(-q * 2 / 3., t2LVec, genVertex_, -q * 6, 3, false);
438  reco::GenParticle w2(-q, w2LVec, genVertex_, -q * 24, 3, true);
439  reco::GenParticle b2(0, bJet2.p4(), genVertex_, -q * 5, 1, true);
440  reco::GenParticle u2(0, wJet1.p4(), genVertex_, -2 * q, 1, true);
441  reco::GenParticle d2(0, wJet2.p4(), genVertex_, q, 1, true);
442 
443  pseudoTop->push_back(t1);
444  pseudoTop->push_back(t2);
445 
446  pseudoTop->push_back(w1);
447  pseudoTop->push_back(b1);
448  pseudoTop->push_back(l1);
449  pseudoTop->push_back(n1);
450 
451  pseudoTop->push_back(w2);
452  pseudoTop->push_back(b2);
453  pseudoTop->push_back(u2);
454  pseudoTop->push_back(d2);
455  }
456  } while (false);
457 
458  if (pseudoTop->size() == 10) // If pseudtop decay tree is completed
459  {
460  // t->W+b
461  pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 2)); // t->W
462  pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 3)); // t->b
463  pseudoTop->at(2).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->W
464  pseudoTop->at(3).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->b
465 
466  // W->lv or W->jj
467  pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 4));
468  pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 5));
469  pseudoTop->at(4).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
470  pseudoTop->at(5).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
471 
472  // tbar->W-b
473  pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 6));
474  pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 7));
475  pseudoTop->at(6).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
476  pseudoTop->at(7).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
477 
478  // W->jj
479  pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 8));
480  pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 9));
481  pseudoTop->at(8).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
482  pseudoTop->at(9).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
483  }
484 
485  event.put(std::move(neutrinos), "neutrinos");
486  event.put(std::move(leptons), "leptons");
487  event.put(std::move(jets), "jets");
488 
489  event.put(std::move(pseudoTop));
490 }

References funct::abs(), b1, b2, RazorAnalyzer::bJet1, RazorAnalyzer::bJet2, reco::Candidate::charge(), funct::cos(), PbPb_ZMuSkimMuonDPG_cff::deltaR, symbols::dm, HGC3DClusterGenMatchSelector_cfi::dR, MillePedeFileConverter_cfg::e, finalStateToken_, fjJetDef_, fjLepDef_, genParticleToken_, genVertex_, mps_fire::i, isBHadron(), isFromHadron(), CommonMethods::isnan(), edm::Ptr< T >::isNonnull(), edm::Ptr< T >::isNull(), dqmiolumiharvest::j, singleTopDQM_cfi::jets, HLT_FULL_cff::leptons, visualization-live-secondInstance_cfg::m, EgHLTOffHistBins_cfi::mass, maxJetEta_, maxLeptonEta_, maxLeptonEtaDilepton_, maxLeptonEtaSemilepton_, maxVetoLeptonEtaSemilepton_, minDileptonMassDilepton_, HLT_FULL_cff::minDR, minJetPt_, minLeptonPt_, minLeptonPtDilepton_, minLeptonPtSemilepton_, minMETSemiLepton_, minMtWSemiLepton_, minVetoLeptonPtSemilepton_, eostools::move(), dqmiodumpmetadata::n, TtSemiLepHypHitFit_cfi::neutrinos, GetRecoTauVFromDQM_MC_cff::next, or, AlCaHLTBitMon_ParallelJobs::p, reco::Candidate::pdgId(), pseudoTop_cfi::pseudoTop, reco::Candidate::pt(), submitPVResolutionJobs::q, q1, q2, L1EGammaCrystalsEmulatorProducer_cfi::scale, jetUpdater_cfi::sort, mathSSE::sqrt(), mps_update::status, RandomServiceHelper::t1, RandomServiceHelper::t2, tMass_, MetAnalyzer::u2, w2, wMass_, and reco::writeSpecific().

Member Data Documentation

◆ finalStateToken_

const edm::EDGetTokenT<edm::View<reco::Candidate> > PseudoTopProducer::finalStateToken_
private

Definition at line 33 of file PseudoTopProducer.h.

Referenced by produce().

◆ fjJetDef_

std::shared_ptr<JetDef> PseudoTopProducer::fjJetDef_
private

Definition at line 44 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

◆ fjLepDef_

std::shared_ptr<JetDef> PseudoTopProducer::fjLepDef_
private

Definition at line 44 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

◆ genParticleToken_

const edm::EDGetTokenT<edm::View<reco::Candidate> > PseudoTopProducer::genParticleToken_
private

Definition at line 34 of file PseudoTopProducer.h.

Referenced by produce().

◆ genVertex_

reco::Particle::Point PseudoTopProducer::genVertex_
private

Definition at line 45 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

◆ maxJetEta_

const double PseudoTopProducer::maxJetEta_
private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

◆ maxLeptonEta_

const double PseudoTopProducer::maxLeptonEta_
private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

◆ maxLeptonEtaDilepton_

const double PseudoTopProducer::maxLeptonEtaDilepton_
private

Definition at line 37 of file PseudoTopProducer.h.

Referenced by produce().

◆ maxLeptonEtaSemilepton_

const double PseudoTopProducer::maxLeptonEtaSemilepton_
private

Definition at line 39 of file PseudoTopProducer.h.

Referenced by produce().

◆ maxVetoLeptonEtaSemilepton_

const double PseudoTopProducer::maxVetoLeptonEtaSemilepton_
private

Definition at line 40 of file PseudoTopProducer.h.

Referenced by produce().

◆ minDileptonMassDilepton_

const double PseudoTopProducer::minDileptonMassDilepton_
private

Definition at line 38 of file PseudoTopProducer.h.

Referenced by produce().

◆ minJetPt_

const double PseudoTopProducer::minJetPt_
private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

◆ minLeptonPt_

const double PseudoTopProducer::minLeptonPt_
private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

◆ minLeptonPtDilepton_

const double PseudoTopProducer::minLeptonPtDilepton_
private

Definition at line 37 of file PseudoTopProducer.h.

Referenced by produce().

◆ minLeptonPtSemilepton_

const double PseudoTopProducer::minLeptonPtSemilepton_
private

Definition at line 39 of file PseudoTopProducer.h.

Referenced by produce().

◆ minMETSemiLepton_

const double PseudoTopProducer::minMETSemiLepton_
private

Definition at line 41 of file PseudoTopProducer.h.

Referenced by produce().

◆ minMtWSemiLepton_

const double PseudoTopProducer::minMtWSemiLepton_
private

Definition at line 41 of file PseudoTopProducer.h.

Referenced by produce().

◆ minVetoLeptonPtSemilepton_

const double PseudoTopProducer::minVetoLeptonPtSemilepton_
private

Definition at line 40 of file PseudoTopProducer.h.

Referenced by produce().

◆ tMass_

const double PseudoTopProducer::tMass_
private

Definition at line 36 of file PseudoTopProducer.h.

Referenced by produce().

◆ wMass_

const double PseudoTopProducer::wMass_
private

Definition at line 36 of file PseudoTopProducer.h.

Referenced by produce().

PseudoTopProducer::fjJetDef_
std::shared_ptr< JetDef > fjJetDef_
Definition: PseudoTopProducer.h:44
RandomServiceHelper.t2
t2
Definition: RandomServiceHelper.py:257
reco::GenParticleRef
edm::Ref< GenParticleCollection > GenParticleRef
persistent reference to a GenParticle
Definition: GenParticleFwd.h:17
w2
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
reco::writeSpecific
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:34
mps_fire.i
i
Definition: mps_fire.py:428
GreaterByPt
Definition: PtComparator.h:24
PseudoTopProducer::genParticleToken_
const edm::EDGetTokenT< edm::View< reco::Candidate > > genParticleToken_
Definition: PseudoTopProducer.h:34
reco::GenJet
Jets made from MC generator particles.
Definition: GenJet.h:23
PseudoTopProducer::minLeptonPtDilepton_
const double minLeptonPtDilepton_
Definition: PseudoTopProducer.h:37
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
edm::RefProd::id
ProductID id() const
Accessor for product ID.
Definition: RefProd.h:124
HLT_FULL_cff.minDR
minDR
Definition: HLT_FULL_cff.py:20381
L1EGammaCrystalsEmulatorProducer_cfi.scale
scale
Definition: L1EGammaCrystalsEmulatorProducer_cfi.py:10
ecalCompactTrigPrim_cfi.outColl
outColl
Definition: ecalCompactTrigPrim_cfi.py:6
PseudoTopProducer::genVertex_
reco::Particle::Point genVertex_
Definition: PseudoTopProducer.h:45
PseudoTopProducer::maxLeptonEtaDilepton_
const double maxLeptonEtaDilepton_
Definition: PseudoTopProducer.h:37
PseudoTopProducer::isBHadron
bool isBHadron(const reco::Candidate *p) const
Definition: PseudoTopProducer.cc:518
reco::GenJetCollection
std::vector< GenJet > GenJetCollection
collection of GenJet objects
Definition: GenJetCollection.h:14
HLT_FULL_cff.leptons
leptons
Definition: HLT_FULL_cff.py:26281
reco::GenParticle
Definition: GenParticle.h:21
mps_update.status
status
Definition: mps_update.py:69
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
PseudoTopProducer::maxLeptonEta_
const double maxLeptonEta_
Definition: PseudoTopProducer.h:35
reco::GenParticleCollection
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
Definition: GenParticleFwd.h:13
reco::Candidate::pt
virtual double pt() const =0
transverse momentum
PseudoTopProducer::maxJetEta_
const double maxJetEta_
Definition: PseudoTopProducer.h:35
reco::Candidate::status
virtual int status() const =0
status word
PseudoTopProducer::maxLeptonEtaSemilepton_
const double maxLeptonEtaSemilepton_
Definition: PseudoTopProducer.h:39
PseudoTopProducer::minMETSemiLepton_
const double minMETSemiLepton_
Definition: PseudoTopProducer.h:41
PseudoTopProducer::minMtWSemiLepton_
const double minMtWSemiLepton_
Definition: PseudoTopProducer.h:41
singleTopDQM_cfi.jets
jets
Definition: singleTopDQM_cfi.py:42
b2
static constexpr float b2
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
PseudoTopProducer::wMass_
const double wMass_
Definition: PseudoTopProducer.h:36
PseudoTopProducer::minVetoLeptonPtSemilepton_
const double minVetoLeptonPtSemilepton_
Definition: PseudoTopProducer.h:40
edm::Handle
Definition: AssociativeIterator.h:50
PseudoTopProducer::tMass_
const double tMass_
Definition: PseudoTopProducer.h:36
edm::Ref< GenParticleCollection >
pseudoTop_cfi.leptonConeSize
leptonConeSize
Definition: pseudoTop_cfi.py:10
PseudoTopProducer::minLeptonPtSemilepton_
const double minLeptonPtSemilepton_
Definition: PseudoTopProducer.h:39
b1
static constexpr float b1
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
reco::Candidate::numberOfMothers
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
reco::Candidate::charge
virtual int charge() const =0
electric charge
pseudoTop_cfi.pseudoTop
pseudoTop
Definition: pseudoTop_cfi.py:3
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:72
PseudoTopProducer::fjLepDef_
std::shared_ptr< JetDef > fjLepDef_
Definition: PseudoTopProducer.h:44
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CommonMethods.isnan
def isnan(num)
Definition: CommonMethods.py:98
q2
double q2[4]
Definition: TauolaWrapper.h:88
PbPb_ZMuSkimMuonDPG_cff.deltaR
deltaR
Definition: PbPb_ZMuSkimMuonDPG_cff.py:63
PseudoTopProducer::getLast
const reco::Candidate * getLast(const reco::Candidate *p)
Definition: PseudoTopProducer.cc:492
edm::View
Definition: CaloClusterFwd.h:14
q1
double q1[4]
Definition: TauolaWrapper.h:87
RazorAnalyzer.bJet1
bJet1
do the razor with one or two b jets (medium CSV)
Definition: RazorAnalyzer.py:234
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
cand
Definition: decayParser.h:32
symbols.dm
dm
Definition: symbols.py:66
PseudoTopProducer::maxVetoLeptonEtaSemilepton_
const double maxVetoLeptonEtaSemilepton_
Definition: PseudoTopProducer.h:40
EgammaValidation_cff.pdgId
pdgId
Definition: EgammaValidation_cff.py:118
submitPVResolutionJobs.q
q
Definition: submitPVResolutionJobs.py:84
PseudoTopProducer::minLeptonPt_
const double minLeptonPt_
Definition: PseudoTopProducer.h:35
reco::Candidate::pdgId
virtual int pdgId() const =0
PDG identifier.
edm::Ptr< Candidate >
reco::Candidate
Definition: Candidate.h:27
RazorAnalyzer.bJet2
bJet2
Definition: RazorAnalyzer.py:235
PseudoTopProducer::isFromHadron
bool isFromHadron(const reco::Candidate *p) const
Definition: PseudoTopProducer.cc:501
TtSemiLepHypHitFit_cfi.neutrinos
neutrinos
Definition: TtSemiLepHypHitFit_cfi.py:17
PseudoTopProducer::finalStateToken_
const edm::EDGetTokenT< edm::View< reco::Candidate > > finalStateToken_
Definition: PseudoTopProducer.h:33
eostools.move
def move(src, dest)
Definition: eostools.py:511
particleLevel_cfi.jetConeSize
jetConeSize
Definition: particleLevel_cfi.py:17
LorentzVector
math::XYZTLorentzVector LorentzVector
Definition: HLTMuonMatchAndPlot.h:49
PseudoTopProducer::minDileptonMassDilepton_
const double minDileptonMassDilepton_
Definition: PseudoTopProducer.h:38
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
or
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
edm::Ptr::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
reco::Particle::Point
math::XYZPoint Point
point in the space
Definition: Particle.h:25
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
PseudoTopProducer::minJetPt_
const double minJetPt_
Definition: PseudoTopProducer.h:35
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
reco::Candidate::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:36
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
MetAnalyzer.u2
u2
Definition: MetAnalyzer.py:61
edm::Ptr::isNull
bool isNull() const
Checks for null.
Definition: Ptr.h:142
edm::InputTag
Definition: InputTag.h:15
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16