CMS 3D CMS Logo

List of all members | Public Member Functions | Private Types | Private Member Functions | Private Attributes
PseudoTopProducer Class Reference
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
 
 EDProducer (const EDProducer &)=delete
 
bool hasAbilityToProduceInBeginLumis () const final
 
bool hasAbilityToProduceInBeginProcessBlocks () const final
 
bool hasAbilityToProduceInBeginRuns () const final
 
bool hasAbilityToProduceInEndLumis () const final
 
bool hasAbilityToProduceInEndProcessBlocks () const final
 
bool hasAbilityToProduceInEndRuns () const final
 
const EDProduceroperator= (const EDProducer &)=delete
 

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::CandidategetLast (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< JetDeffjJetDef_
 
std::shared_ptr< JetDeffjLepDef_
 
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<>
using CacheTypes = CacheContexts< T... >
 
using GlobalCache = typename CacheTypes::GlobalCache
 
using HasAbility = AbilityChecker< T... >
 
using InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache
 
using LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache
 
using LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >
 
using LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache
 
using RunCache = typename CacheTypes::RunCache
 
using RunContext = RunContextT< RunCache, GlobalCache >
 
using RunSummaryCache = typename CacheTypes::RunSummaryCache
 

Detailed Description

Definition at line 18 of file PseudoTopProducer.cc.

Member Typedef Documentation

◆ JetDef

typedef fastjet::JetDefinition PseudoTopProducer::JetDef
private

Definition at line 46 of file PseudoTopProducer.cc.

◆ LorentzVector

Definition at line 33 of file PseudoTopProducer.cc.

Constructor & Destructor Documentation

◆ PseudoTopProducer()

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

Definition at line 55 of file PseudoTopProducer.cc.

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

56  : finalStateToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("finalStates"))),
57  genParticleToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("genParticles"))),
58  minLeptonPt_(pset.getParameter<double>("minLeptonPt")),
59  maxLeptonEta_(pset.getParameter<double>("maxLeptonEta")),
60  minJetPt_(pset.getParameter<double>("minJetPt")),
61  maxJetEta_(pset.getParameter<double>("maxJetEta")),
62  wMass_(pset.getParameter<double>("wMass")),
63  tMass_(pset.getParameter<double>("tMass")),
64  minLeptonPtDilepton_(pset.getParameter<double>("minLeptonPtDilepton")),
65  maxLeptonEtaDilepton_(pset.getParameter<double>("maxLeptonEtaDilepton")),
66  minDileptonMassDilepton_(pset.getParameter<double>("minDileptonMassDilepton")),
67  minLeptonPtSemilepton_(pset.getParameter<double>("minLeptonPtSemilepton")),
68  maxLeptonEtaSemilepton_(pset.getParameter<double>("maxLeptonEtaSemilepton")),
69  minVetoLeptonPtSemilepton_(pset.getParameter<double>("minVetoLeptonPtSemilepton")),
70  maxVetoLeptonEtaSemilepton_(pset.getParameter<double>("maxVetoLeptonEtaSemilepton")),
71  minMETSemiLepton_(pset.getParameter<double>("minMETSemiLepton")),
72  minMtWSemiLepton_(pset.getParameter<double>("minMtWSemiLepton")) {
73  const double leptonConeSize = pset.getParameter<double>("leptonConeSize");
74  const double jetConeSize = pset.getParameter<double>("jetConeSize");
75  fjLepDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, leptonConeSize);
76  fjJetDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, jetConeSize);
77 
79 
80  produces<reco::GenParticleCollection>("neutrinos");
81  produces<reco::GenJetCollection>("leptons");
82  produces<reco::GenJetCollection>("jets");
83 
84  produces<reco::GenParticleCollection>();
85 }
const edm::EDGetTokenT< edm::View< reco::Candidate > > genParticleToken_
reco::Particle::Point genVertex_
const double maxLeptonEtaSemilepton_
const double maxLeptonEta_
const double minVetoLeptonPtSemilepton_
const double maxLeptonEtaDilepton_
std::shared_ptr< JetDef > fjJetDef_
std::shared_ptr< JetDef > fjLepDef_
const double minMtWSemiLepton_
const double minMETSemiLepton_
const double maxJetEta_
const double minLeptonPtSemilepton_
math::XYZPoint Point
point in the space
Definition: Particle.h:25
const double maxVetoLeptonEtaSemilepton_
const double minLeptonPt_
const double minDileptonMassDilepton_
const edm::EDGetTokenT< edm::View< reco::Candidate > > finalStateToken_
const double minLeptonPtDilepton_
const double minJetPt_

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 599 of file PseudoTopProducer.cc.

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().

601  {
602  reco::GenParticle pOut(*dynamic_cast<const reco::GenParticle*>(p));
603  pOut.clearMothers();
604  pOut.clearDaughters();
605  pOut.resetMothers(refHandle.id());
606  pOut.resetDaughters(refHandle.id());
607 
608  outColl->push_back(pOut);
609 
610  return reco::GenParticleRef(refHandle, outColl->size() - 1);
611 }
edm::Ref< GenParticleCollection > GenParticleRef
persistent reference to a GenParticle
ProductID id() const
Accessor for product ID.
Definition: RefProd.h:124

◆ getLast()

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

Definition at line 534 of file PseudoTopProducer.cc.

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

534  {
535  for (size_t i = 0, n = p->numberOfDaughters(); i < n; ++i) {
536  const reco::Candidate* dau = p->daughter(i);
537  if (p->pdgId() == dau->pdgId())
538  return getLast(dau);
539  }
540  return p;
541 }
const reco::Candidate * getLast(const reco::Candidate *p)
virtual int pdgId() const =0
PDG identifier.

◆ 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 560 of file PseudoTopProducer.cc.

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

Referenced by produce().

560  {
561  const unsigned int absPdgId = abs(p->pdgId());
562  if (!isBHadron(absPdgId))
563  return false;
564 
565  // Do not consider this particle if it has B hadron daughter
566  // For example, B* -> B0 + photon; then we drop B* and take B0 only
567  for (int i = 0, n = p->numberOfDaughters(); i < n; ++i) {
568  const reco::Candidate* dau = p->daughter(i);
569  if (isBHadron(abs(dau->pdgId())))
570  return false;
571  }
572 
573  return true;
574 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual int pdgId() const =0
PDG identifier.
bool isBHadron(const reco::Candidate *p) const

◆ isBHadron() [2/2]

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

Definition at line 576 of file PseudoTopProducer.cc.

576  {
577  if (absPdgId <= 100)
578  return false; // Fundamental particles and MC internals
579  if (absPdgId >= 1000000000)
580  return false; // Nuclei, +-10LZZZAAAI
581 
582  // General form of PDG ID is 7 digit form
583  // +- n nr nL nq1 nq2 nq3 nJ
584  //const int nJ = absPdgId % 10; // Spin
585  const int nq3 = (absPdgId / 10) % 10;
586  const int nq2 = (absPdgId / 100) % 10;
587  const int nq1 = (absPdgId / 1000) % 10;
588 
589  if (nq3 == 0)
590  return false; // Diquarks
591  if (nq1 == 0 and nq2 == 5)
592  return true; // B mesons
593  if (nq1 == 5)
594  return true; // B baryons
595 
596  return false;
597 }

◆ isFromHadron()

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

Definition at line 543 of file PseudoTopProducer.cc.

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().

543  {
544  for (size_t i = 0, n = p->numberOfMothers(); i < n; ++i) {
545  const reco::Candidate* mother = p->mother(i);
546  if (!mother)
547  continue;
548  if (mother->status() == 4 or mother->numberOfMothers() == 0)
549  continue; // Skip incident beam
550  const int pdgId = abs(mother->pdgId());
551 
552  if (pdgId > 100)
553  return true;
554  else if (isFromHadron(mother))
555  return true;
556  }
557  return false;
558 }
virtual int status() const =0
status word
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
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
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual int pdgId() const =0
PDG identifier.
bool isFromHadron(const reco::Candidate *p) const

◆ produce()

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

Build jets

Definition at line 87 of file PseudoTopProducer.cc.

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, PDWG_EXODelayedJetMET_cff::jets, HLT_2022v15_cff::leptons, visualization-live-secondInstance_cfg::m, EgHLTOffHistBins_cfi::mass, maxJetEta_, maxLeptonEta_, maxLeptonEtaDilepton_, maxLeptonEtaSemilepton_, maxVetoLeptonEtaSemilepton_, minDileptonMassDilepton_, HLT_2022v15_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, L1EGammaClusterEmuProducer_cfi::scale, jetsAK4_CHS_cff::sort, mathSSE::sqrt(), mps_update::status, RandomServiceHelper::t1, RandomServiceHelper::t2, tMass_, MetAnalyzer::u2, w1, w2, wMass_, and reco::writeSpecific().

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

Member Data Documentation

◆ finalStateToken_

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

Definition at line 36 of file PseudoTopProducer.cc.

Referenced by produce().

◆ fjJetDef_

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

Definition at line 47 of file PseudoTopProducer.cc.

Referenced by produce(), and PseudoTopProducer().

◆ fjLepDef_

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

Definition at line 47 of file PseudoTopProducer.cc.

Referenced by produce(), and PseudoTopProducer().

◆ genParticleToken_

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

Definition at line 37 of file PseudoTopProducer.cc.

Referenced by produce().

◆ genVertex_

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

Definition at line 48 of file PseudoTopProducer.cc.

Referenced by produce(), and PseudoTopProducer().

◆ maxJetEta_

const double PseudoTopProducer::maxJetEta_
private

Definition at line 38 of file PseudoTopProducer.cc.

Referenced by produce().

◆ maxLeptonEta_

const double PseudoTopProducer::maxLeptonEta_
private

Definition at line 38 of file PseudoTopProducer.cc.

Referenced by produce().

◆ maxLeptonEtaDilepton_

const double PseudoTopProducer::maxLeptonEtaDilepton_
private

Definition at line 40 of file PseudoTopProducer.cc.

Referenced by produce().

◆ maxLeptonEtaSemilepton_

const double PseudoTopProducer::maxLeptonEtaSemilepton_
private

Definition at line 42 of file PseudoTopProducer.cc.

Referenced by produce().

◆ maxVetoLeptonEtaSemilepton_

const double PseudoTopProducer::maxVetoLeptonEtaSemilepton_
private

Definition at line 43 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minDileptonMassDilepton_

const double PseudoTopProducer::minDileptonMassDilepton_
private

Definition at line 41 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minJetPt_

const double PseudoTopProducer::minJetPt_
private

Definition at line 38 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minLeptonPt_

const double PseudoTopProducer::minLeptonPt_
private

Definition at line 38 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minLeptonPtDilepton_

const double PseudoTopProducer::minLeptonPtDilepton_
private

Definition at line 40 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minLeptonPtSemilepton_

const double PseudoTopProducer::minLeptonPtSemilepton_
private

Definition at line 42 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minMETSemiLepton_

const double PseudoTopProducer::minMETSemiLepton_
private

Definition at line 44 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minMtWSemiLepton_

const double PseudoTopProducer::minMtWSemiLepton_
private

Definition at line 44 of file PseudoTopProducer.cc.

Referenced by produce().

◆ minVetoLeptonPtSemilepton_

const double PseudoTopProducer::minVetoLeptonPtSemilepton_
private

Definition at line 43 of file PseudoTopProducer.cc.

Referenced by produce().

◆ tMass_

const double PseudoTopProducer::tMass_
private

Definition at line 39 of file PseudoTopProducer.cc.

Referenced by produce().

◆ wMass_

const double PseudoTopProducer::wMass_
private

Definition at line 39 of file PseudoTopProducer.cc.

Referenced by produce().