CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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::EDProducer edm::ProducerBase edm::EDConsumerBase edm::ProductRegistryHelper

Public Member Functions

void produce (edm::Event &event, const edm::EventSetup &eventSetup) override
 
 PseudoTopProducer (const edm::ParameterSet &pset)
 
- Public Member Functions inherited from edm::EDProducer
 EDProducer ()
 
ModuleDescription const & moduleDescription () const
 
virtual ~EDProducer ()
 
- Public Member Functions inherited from edm::ProducerBase
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 ProducerBase ()
 
void registerProducts (ProducerBase *, ProductRegistry *, ModuleDescription const &)
 
std::function< void(BranchDescription
const &)> 
registrationCallback () const
 used by the fwk to register list of products More...
 
virtual ~ProducerBase ()
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfoconsumesInfo () const
 
 EDConsumerBase ()
 
ProductHolderIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
std::vector
< ProductHolderIndexAndSkipBit >
const & 
itemsToGetFromEvent () const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesDependentUpon (std::string const &iProcessName, std::string const &iModuleLabel, bool iPrint, std::vector< char const * > &oModuleLabels) const
 
void modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
 
bool registeredToConsume (ProductHolderIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
 
virtual ~EDConsumerBase ()
 

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)
 
bool isBHadron (const reco::Candidate *p) const
 
bool isBHadron (const unsigned int pdgId) const
 
bool isFromHadron (const reco::Candidate *p) const
 

Private Attributes

edm::EDGetTokenT< edm::View
< reco::Candidate > > 
finalStateToken_
 
std::shared_ptr< JetDeffjJetDef_
 
std::shared_ptr< JetDeffjLepDef_
 
edm::EDGetTokenT< edm::View
< reco::Candidate > > 
genParticleToken_
 
reco::Particle::Point genVertex_
 
const double jetMaxEta_
 
const double jetMinPt_
 
const double leptonMaxEta_
 
const double leptonMinPt_
 
const double tMass_
 
const double wMass_
 

Additional Inherited Members

- Public Types inherited from edm::EDProducer
typedef EDProducer ModuleType
 
- Public Types inherited from edm::ProducerBase
typedef
ProductRegistryHelper::TypeLabelList 
TypeLabelList
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::EDProducer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &descriptions)
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 

Detailed Description

Definition at line 13 of file PseudoTopProducer.h.

Member Typedef Documentation

typedef fastjet::JetDefinition PseudoTopProducer::JetDef
private

Definition at line 36 of file PseudoTopProducer.h.

Definition at line 28 of file PseudoTopProducer.h.

Constructor & Destructor Documentation

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

Definition at line 12 of file PseudoTopProducer.cc.

References finalStateToken_, fjJetDef_, fjLepDef_, genParticleToken_, genVertex_, and edm::ParameterSet::getParameter().

12  :
13  leptonMinPt_(pset.getParameter<double>("leptonMinPt")),
14  leptonMaxEta_(pset.getParameter<double>("leptonMaxEta")),
15  jetMinPt_(pset.getParameter<double>("jetMinPt")),
16  jetMaxEta_(pset.getParameter<double>("jetMaxEta")),
17  wMass_(pset.getParameter<double>("wMass")),
18  tMass_(pset.getParameter<double>("tMass"))
19 {
20  finalStateToken_ = consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("finalStates"));
21  genParticleToken_ = consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("genParticles"));
22 
23  const double leptonConeSize = pset.getParameter<double>("leptonConeSize");
24  const double jetConeSize = pset.getParameter<double>("jetConeSize");
25  fjLepDef_ = std::shared_ptr<JetDef>(new JetDef(fastjet::antikt_algorithm, leptonConeSize));
26  fjJetDef_ = std::shared_ptr<JetDef>(new JetDef(fastjet::antikt_algorithm, jetConeSize));
27 
29 
30  produces<reco::GenParticleCollection>("neutrinos");
31  produces<reco::GenJetCollection>("leptons");
32  produces<reco::GenJetCollection>("jets");
33 
34  produces<reco::GenParticleCollection>();
35 
36 }
const double leptonMaxEta_
T getParameter(std::string const &) const
reco::Particle::Point genVertex_
const double leptonMinPt_
const double jetMaxEta_
edm::EDGetTokenT< edm::View< reco::Candidate > > finalStateToken_
fastjet::JetDefinition JetDef
std::shared_ptr< JetDef > fjJetDef_
std::shared_ptr< JetDef > fjLepDef_
math::XYZPoint Point
point in the space
Definition: Particle.h:25
edm::EDGetTokenT< edm::View< reco::Candidate > > genParticleToken_
const double jetMinPt_

Member Function Documentation

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

Definition at line 563 of file PseudoTopProducer.cc.

References reco::CompositeRefCandidateT< D >::clearDaughters(), reco::CompositeRefCandidateT< D >::clearMothers(), edm::RefProd< T >::id(), reco::CompositeRefCandidateT< D >::resetDaughters(), and reco::CompositeRefCandidateT< D >::resetMothers().

565 {
566  reco::GenParticle pOut(*dynamic_cast<const reco::GenParticle*>(p));
567  pOut.clearMothers();
568  pOut.clearDaughters();
569  pOut.resetMothers(refHandle.id());
570  pOut.resetDaughters(refHandle.id());
571 
572  outColl->push_back(pOut);
573 
574  return reco::GenParticleRef(refHandle, outColl->size()-1);
575 }
edm::Ref< GenParticleCollection > GenParticleRef
persistent reference to a GenParticle
ProductID id() const
Accessor for product ID.
Definition: RefProd.h:137
const reco::Candidate * PseudoTopProducer::getLast ( const reco::Candidate p)
private

Definition at line 504 of file PseudoTopProducer.cc.

References reco::Candidate::daughter(), i, gen::n, reco::Candidate::numberOfDaughters(), AlCaHLTBitMon_ParallelJobs::p, and reco::Candidate::pdgId().

505 {
506  for ( size_t i=0, n=p->numberOfDaughters(); i<n; ++i )
507  {
508  const reco::Candidate* dau = p->daughter(i);
509  if ( p->pdgId() == dau->pdgId() ) return getLast(dau);
510  }
511  return p;
512 }
int i
Definition: DBlmapReader.cc:9
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
virtual size_type numberOfDaughters() const =0
number of daughters
const reco::Candidate * getLast(const reco::Candidate *p)
virtual int pdgId() const =0
PDG identifier.
bool PseudoTopProducer::isBHadron ( const reco::Candidate p) const
private

Definition at line 528 of file PseudoTopProducer.cc.

References funct::abs(), reco::Candidate::daughter(), i, gen::n, reco::Candidate::numberOfDaughters(), and reco::Candidate::pdgId().

Referenced by produce().

529 {
530  const unsigned int absPdgId = abs(p->pdgId());
531  if ( !isBHadron(absPdgId) ) return false;
532 
533  // Do not consider this particle if it has B hadron daughter
534  // For example, B* -> B0 + photon; then we drop B* and take B0 only
535  for ( int i=0, n=p->numberOfDaughters(); i<n; ++i )
536  {
537  const reco::Candidate* dau = p->daughter(i);
538  if ( isBHadron(abs(dau->pdgId())) ) return false;
539  }
540 
541  return true;
542 }
bool isBHadron(const reco::Candidate *p) const
int i
Definition: DBlmapReader.cc:9
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
virtual size_type numberOfDaughters() const =0
number of daughters
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual int pdgId() const =0
PDG identifier.
bool PseudoTopProducer::isBHadron ( const unsigned int  pdgId) const
private

Definition at line 544 of file PseudoTopProducer.cc.

545 {
546  if ( absPdgId <= 100 ) return false; // Fundamental particles and MC internals
547  if ( absPdgId >= 1000000000 ) return false; // Nuclei, +-10LZZZAAAI
548 
549  // General form of PDG ID is 7 digit form
550  // +- n nr nL nq1 nq2 nq3 nJ
551  //const int nJ = absPdgId % 10; // Spin
552  const int nq3 = (absPdgId / 10) % 10;
553  const int nq2 = (absPdgId / 100) % 10;
554  const int nq1 = (absPdgId / 1000) % 10;
555 
556  if ( nq3 == 0 ) return false; // Diquarks
557  if ( nq1 == 0 and nq2 == 5 ) return true; // B mesons
558  if ( nq1 == 5 ) return true; // B baryons
559 
560  return false;
561 }
bool PseudoTopProducer::isFromHadron ( const reco::Candidate p) const
private

Definition at line 514 of file PseudoTopProducer.cc.

References funct::abs(), i, reco::Candidate::mother(), gen::n, reco::Candidate::numberOfMothers(), SingleMuPt40Fwdv2_cfi_GEN_SIM::pdgId, and reco::Candidate::pdgId().

Referenced by produce().

515 {
516  for ( size_t i=0, n=p->numberOfMothers(); i<n; ++i )
517  {
518  const reco::Candidate* mother = p->mother(i);
519  if ( mother->numberOfMothers() == 0 ) continue; // Skip incident beam
520  const int pdgId = abs(mother->pdgId());
521 
522  if ( pdgId > 100 ) return true;
523  else if ( isFromHadron(mother) ) return true;
524  }
525  return false;
526 }
int i
Definition: DBlmapReader.cc:9
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
bool isFromHadron(const reco::Candidate *p) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual int pdgId() const =0
PDG identifier.
void PseudoTopProducer::produce ( edm::Event event,
const edm::EventSetup eventSetup 
)
overridevirtual

Build jets

Implements edm::EDProducer.

Definition at line 38 of file PseudoTopProducer.cc.

References funct::abs(), reco::CompositeRefCandidateT< D >::addDaughter(), RazorAnalyzer::bJet1, RazorAnalyzer::bJet2, reco::Candidate::charge(), symbols::dm, alignCSCRings::e, reco::Candidate::energy(), finalStateToken_, fjJetDef_, fjLepDef_, genParticleToken_, genVertex_, i, cmsHarvester::index, isBHadron(), isFromHadron(), edm::detail::isnan(), edm::Ptr< T >::isNonnull(), edm::Ptr< T >::isNull(), j, jetMaxEta_, jetMinPt_, fwrapper::jets, leptonMaxEta_, leptonMinPt_, HLT_FULL_cff::leptons, visualization-live-secondInstance_cfg::m, reco::Candidate::mother(), gen::n, reco::Candidate::numberOfMothers(), or, reco::Candidate::p(), AlCaHLTBitMon_ParallelJobs::p, reco::Candidate::p4(), reco::Candidate::pdgId(), pseudoTop_cfi::pseudoTop, reco::Candidate::pt(), reco::Candidate::px(), reco::Candidate::py(), reco::Candidate::pz(), pileupReCalc_HLTpaths::scale, mps_update::status, reco::Candidate::status(), reco::t2, tMass_, reco::Candidate::vertex(), w2, wMass_, and reco::writeSpecific().

39 {
40  edm::Handle<edm::View<reco::Candidate> > finalStateHandle;
41  event.getByToken(finalStateToken_, finalStateHandle);
42 
43  edm::Handle<edm::View<reco::Candidate> > genParticleHandle;
44  event.getByToken(genParticleToken_, genParticleHandle);
45 
46  std::auto_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
47  std::auto_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
48  std::auto_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
49  auto neutrinosRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("neutrinos");
50  auto leptonsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("leptons");
51  auto jetsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("jets");
52 
53  std::auto_ptr<reco::GenParticleCollection> pseudoTop(new reco::GenParticleCollection);
54  auto pseudoTopRefHandle = event.getRefBeforePut<reco::GenParticleCollection>();
55 
56  // Collect unstable B-hadrons
57  std::set<size_t> bHadronIdxs;
58  for ( size_t i=0, n=genParticleHandle->size(); i<n; ++i )
59  {
60  const reco::Candidate& p = genParticleHandle->at(i);
61  const int status = p.status();
62  if ( status == 1 ) continue;
63 
64  // Collect B-hadrons, to be used in b tagging
65  if ( isBHadron(&p) ) bHadronIdxs.insert(i);
66  }
67 
68  // Collect stable leptons and neutrinos
69  size_t nStables = 0;
70  std::vector<size_t> leptonIdxs;
71  std::set<size_t> neutrinoIdxs;
72  for ( size_t i=0, n=finalStateHandle->size(); i<n; ++i )
73  {
74  const reco::Candidate& p = finalStateHandle->at(i);
75  const int absPdgId = abs(p.pdgId());
76  if ( p.status() != 1 ) continue;
77 
78  ++nStables;
79  if ( p.numberOfMothers() == 0 ) continue; // Skip orphans (if exists)
80  if ( p.mother()->status() == 4 ) continue; // Treat particle as hadronic if directly from the incident beam (protect orphans in MINIAOD)
81  if ( isFromHadron(&p) ) continue;
82  switch ( absPdgId )
83  {
84  case 11: case 13: // Leptons
85  case 22: // Photons
86  leptonIdxs.push_back(i);
87  break;
88  case 12: case 14: case 16:
89  neutrinoIdxs.insert(i);
90  neutrinos->push_back(reco::GenParticle(p.charge(), p.p4(), p.vertex(), p.pdgId(), p.status(), true));
91  break;
92  }
93  }
94 
95  // Sort neutrinos by pT.
96  std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
97 
98  // Make dressed leptons with anti-kt(0.1) algorithm
100  std::vector<fastjet::PseudoJet> fjLepInputs;
101  fjLepInputs.reserve(leptonIdxs.size());
102  for ( auto index : leptonIdxs )
103  {
104  const reco::Candidate& p = finalStateHandle->at(index);
105  if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
106 
107  fjLepInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
108  fjLepInputs.back().set_user_index(index);
109  }
110 
112  fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *fjLepDef_);
113  std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(leptonMinPt_));
114 
116  leptons->reserve(fjLepJets.size());
117  std::set<size_t> lepDauIdxs; // keep lepton constituents to remove from GenJet construction
118  for ( auto& fjJet : fjLepJets )
119  {
120  if ( abs(fjJet.eta()) > leptonMaxEta_ ) continue;
121 
122  // Get jet constituents from fastJet
123  const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
124  // Convert to CandidatePtr
125  std::vector<reco::CandidatePtr> constituents;
126  reco::CandidatePtr lepCand;
127  for ( auto& fjConstituent : fjConstituents )
128  {
129  const size_t index = fjConstituent.user_index();
130  reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
131  const int absPdgId = abs(cand->pdgId());
132  if ( absPdgId == 11 or absPdgId == 13 )
133  {
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  if ( lepCand->pt() < fjJet.pt()/2 ) continue; // Central lepton must be the major component
141 
142  const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
143  reco::GenJet lepJet;
144  reco::writeSpecific(lepJet, jetP4, genVertex_, constituents, eventSetup);
145 
146  lepJet.setPdgId(lepCand->pdgId());
147  lepJet.setCharge(lepCand->charge());
148 
149  const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
150  lepJet.setJetArea(jetArea);
151 
152  leptons->push_back(lepJet);
153 
154  // Keep constituent indices to be used in the next step.
155  for ( auto& fjConstituent : fjConstituents )
156  {
157  lepDauIdxs.insert(fjConstituent.user_index());
158  }
159  }
160 
161  // Now proceed to jets.
162  // Jets: anti-kt excluding the e, mu, nu, and photons in selected leptons.
164  std::vector<fastjet::PseudoJet> fjJetInputs;
165  fjJetInputs.reserve(nStables);
166  for ( size_t i=0, n=finalStateHandle->size(); i<n; ++i )
167  {
168  const reco::Candidate& p = finalStateHandle->at(i);
169  if ( p.status() != 1 ) continue;
170  if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
171 
172  if ( neutrinoIdxs.find(i) != neutrinoIdxs.end() ) continue;
173  if ( lepDauIdxs.find(i) != lepDauIdxs.end() ) continue;
174 
175  fjJetInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
176  fjJetInputs.back().set_user_index(i);
177  }
179  for ( auto index : bHadronIdxs )
180  {
181  const reco::Candidate& p = genParticleHandle->at(index);
182  if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
183 
184  const double scale = 1e-20/p.p();
185  fjJetInputs.push_back(fastjet::PseudoJet(p.px()*scale, p.py()*scale, p.pz()*scale, p.energy()*scale));
186  fjJetInputs.back().set_user_index(index);
187  }
188 
190  fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *fjJetDef_);
191  std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(jetMinPt_));
192 
194  jets->reserve(fjJets.size());
195  std::vector<size_t> bjetIdxs, ljetIdxs;
196  for ( auto& fjJet : fjJets )
197  {
198  if ( abs(fjJet.eta()) > jetMaxEta_ ) continue;
199 
200  // Get jet constituents from fastJet
201  const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
202  // Convert to CandidatePtr
203  std::vector<reco::CandidatePtr> constituents;
204  bool hasBHadron = false;
205  for ( size_t j=0, m=fjConstituents.size(); j<m; ++j )
206  {
207  const size_t index = fjConstituents[j].user_index();
208  if ( bHadronIdxs.find(index) != bHadronIdxs.end() ) hasBHadron = true;
209  reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
210  constituents.push_back(cand);
211  }
212 
213  const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
214  reco::GenJet genJet;
215  reco::writeSpecific(genJet, jetP4, genVertex_, constituents, eventSetup);
216 
217  const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
218  genJet.setJetArea(jetArea);
219  if ( hasBHadron )
220  {
221  genJet.setPdgId(5);
222  bjetIdxs.push_back(jets->size());
223  }
224  else
225  {
226  ljetIdxs.push_back(jets->size());
227  }
228 
229  jets->push_back(genJet);
230  }
231 
232  // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
233  // NOTE : A C++ trick, use do-while instead of long-nested if-statements.
234  do
235  {
236  // Note : we will do dilepton or semilepton channel only
237  const size_t nLepton = leptons->size();
238 
239  // Collect leptonic-decaying W's
240  std::vector<std::pair<int, int> > wCandIdxs;
241  for ( auto lep = leptons->begin(); lep != leptons->end(); ++lep )
242  {
243  const size_t iLep = lep-leptons->begin();
244  for ( auto nu = neutrinos->begin(); nu != neutrinos->end(); ++nu )
245  {
246  //if ( abs(lep->pdgId()+nu->pdgId()) != 1 ) continue; // Enforce to conserve flavour, this reduces combinatorial bkg
247  const double m = (lep->p4()+nu->p4()).mass();
248  if ( m > 300 ) continue; // Raw mass cut, reduce combinatorial background
249 
250  const size_t iNu = nu-neutrinos->begin();
251  wCandIdxs.push_back(make_pair(iLep, iNu));
252  }
253  }
254 
255  if ( nLepton == 0 or wCandIdxs.empty() ) break; // Skip full hadronic channel
256  else if ( nLepton == 1 and wCandIdxs.size() >= 1 ) // Semi-leptonic channel
257  {
258  double dm = 1e9; // Note: this will be re-used later.
259  int bestLepIdx = -1, bestNuIdx = -1;
260  for ( auto wCandIdx : wCandIdxs )
261  {
262  const int lepIdx = wCandIdx.first;
263  const int nuIdx = wCandIdx.second;
264  const LorentzVector& lepLVec = leptons->at(lepIdx).p4();
265  const LorentzVector& nuLVec = neutrinos->at(nuIdx).p4();
266  const double mW = (lepLVec + nuLVec).mass();
267  if ( mW > 300 ) continue;
268 
269  const double dmNew = abs(mW-wMass_);
270  if ( dmNew < dm )
271  {
272  dm = dmNew;
273  bestLepIdx = lepIdx;
274  bestNuIdx = nuIdx;
275  }
276  }
277  if ( bestLepIdx < 0 ) break; // Actually this should never happen
278  const auto& lepton = leptons->at(bestLepIdx);
279  const auto& neutrino = neutrinos->at(bestNuIdx);
280  const LorentzVector w1LVec = lepton.p4()+neutrino.p4();
281 
282  // Continue to hadronic W
283  // We also collect b jets to be used in the next step.
284  dm = 1e9; // Reset dM to very large value.
285  int bestJ1Idx = -1, bestJ2Idx = -1;
286  for ( auto j1Idx = ljetIdxs.begin(); j1Idx != ljetIdxs.end(); ++j1Idx )
287  {
288  const auto& j1 = jets->at(*j1Idx);
289  for ( auto j2Idx = j1Idx+1; j2Idx != ljetIdxs.end(); ++j2Idx )
290  {
291  const auto& j2 = jets->at(*j2Idx);
292  const double mW = (j1.p4()+j2.p4()).mass();
293  if ( mW > 300 ) continue;
294 
295  const double dmNew = abs(mW-wMass_);
296  if ( dmNew < dm )
297  {
298  dm = dmNew;
299  bestJ1Idx = *j1Idx;
300  bestJ2Idx = *j2Idx;
301  }
302  }
303  }
304  if ( bestJ1Idx < 0 ) break;
305  const auto& wJet1 = jets->at(bestJ1Idx);
306  const auto& wJet2 = jets->at(bestJ2Idx);
307  const LorentzVector w2LVec = wJet1.p4() + wJet2.p4();
308 
309  // Now we have leptonic W and hadronic W.
310  // Contiue to top quarks
311  dm = 1e9; // Reset once again for top combination.
312  int bestB1Idx = -1, bestB2Idx = -1;
313  if ( bjetIdxs.size() < 2 ) break;
314  for ( auto b1Idx : bjetIdxs )
315  {
316  const double t1Mass = (w1LVec + jets->at(b1Idx).p4()).mass();
317  if ( t1Mass > 300 ) continue;
318  for ( auto b2Idx : bjetIdxs )
319  {
320  if ( b1Idx == b2Idx ) continue;
321  const double t2Mass = (w2LVec + jets->at(b2Idx).p4()).mass();
322  if ( t2Mass > 300 ) continue;
323 
324  const double dmNew = abs(t1Mass-tMass_) + abs(t2Mass-tMass_);
325  if ( dmNew < dm )
326  {
327  dm = dmNew;
328  bestB1Idx = b1Idx;
329  bestB2Idx = b2Idx;
330  }
331  }
332  }
333  if ( bestB1Idx < 0 ) break;
334  const auto& bJet1 = jets->at(bestB1Idx);
335  const auto& bJet2 = jets->at(bestB2Idx);
336  const LorentzVector t1LVec = w1LVec + bJet1.p4();
337  const LorentzVector t2LVec = w2LVec + bJet2.p4();
338 
339  // Put all of them into candidate collection
340  if ( true ) // Trick to restrict variables' scope to avoid collision
341  {
342  const int lepQ = lepton.charge();
343 
344  reco::GenParticle t1(lepQ*2/3, t1LVec, genVertex_, lepQ*6, 3, false);
345  reco::GenParticle w1(lepQ, w1LVec, genVertex_, lepQ*24, 3, true);
346  reco::GenParticle b1(0, bJet1.p4(), genVertex_, lepQ*5, 1, true);
347  reco::GenParticle l1(lepQ, lepton.p4(), genVertex_, lepton.pdgId(), 1, true);
348  reco::GenParticle n1(0, neutrino.p4(), genVertex_, neutrino.pdgId(), 1, true);
349 
350  reco::GenParticle t2(-lepQ*2/3, t2LVec, genVertex_, -lepQ*6, 3, false);
351  reco::GenParticle w2(-lepQ, w2LVec, genVertex_, -lepQ*24, 3, true);
352  reco::GenParticle b2(0, bJet2.p4(), genVertex_, -lepQ*5, 1, true);
353  reco::GenParticle qA(0, wJet1.p4(), genVertex_, 1, 1, true);
354  reco::GenParticle qB(0, wJet2.p4(), genVertex_, 2, 1, true);
355 
356  pseudoTop->push_back(t1);
357  pseudoTop->push_back(w1);
358  pseudoTop->push_back(b1);
359  pseudoTop->push_back(l1);
360  pseudoTop->push_back(n1);
361 
362  pseudoTop->push_back(t2);
363  pseudoTop->push_back(w2);
364  pseudoTop->push_back(b2);
365  pseudoTop->push_back(qA);
366  pseudoTop->push_back(qB);
367  }
368  }
369  else if ( nLepton == 2 and wCandIdxs.size() >= 2 ) // Dilepton channel
370  {
371  double dm = 1e9;
372  int bestLep1Idx = -1, bestLep2Idx = -1, bestNu1Idx = -1, bestNu2Idx = -1;
373  for ( auto i : wCandIdxs )
374  {
375  const auto& lepton1 = leptons->at(i.first);
376  if ( lepton1.charge() < 0 ) continue;
377  const auto& neutrino1 = neutrinos->at(i.second);
378  const double mW1 = (lepton1.p4()+neutrino1.p4()).mass();
379  if ( mW1 > 300 ) continue;
380  for ( auto j : wCandIdxs )
381  {
382  if ( i == j ) continue;
383  const auto& lepton2 = leptons->at(j.first);
384  if ( lepton2.charge() > 0 ) continue;
385  const auto& neutrino2 = neutrinos->at(j.second);
386  const double mW2 = (lepton2.p4()+neutrino2.p4()).mass();
387  if ( mW2 > 300 ) continue;
388 
389  const double dmNew = abs(mW1-wMass_) + abs(mW2-wMass_);
390  if ( dmNew < dm )
391  {
392  dm = dmNew;
393  bestLep1Idx = i.first;
394  bestLep2Idx = j.first;
395  bestNu1Idx = i.second;
396  bestNu2Idx = j.second;
397  }
398  }
399  }
400  if ( bestLep1Idx < 0 ) break;
401  const auto& lepton1 = leptons->at(bestLep1Idx);
402  const auto& lepton2 = leptons->at(bestLep2Idx);
403  const auto& neutrino1 = neutrinos->at(bestNu1Idx);
404  const auto& neutrino2 = neutrinos->at(bestNu2Idx);
405  const LorentzVector w1LVec = lepton1.p4()+neutrino1.p4();
406  const LorentzVector w2LVec = lepton2.p4()+neutrino2.p4();
407 
408  // Contiue to top quarks
409  dm = 1e9; // Reset once again for top combination.
410  int bestB1Idx = -1, bestB2Idx = -1;
411  if ( bjetIdxs.size() < 2 ) break;
412  for ( auto b1Idx : bjetIdxs )
413  {
414  const double t1Mass = (w1LVec + jets->at(b1Idx).p4()).mass();
415  if ( t1Mass > 300 ) continue;
416  for ( auto b2Idx : bjetIdxs )
417  {
418  if ( b1Idx == b2Idx ) continue;
419  const double t2Mass = (w2LVec + jets->at(b2Idx).p4()).mass();
420  if ( t2Mass > 300 ) continue;
421 
422  const double dmNew = abs(t1Mass-tMass_) + abs(t2Mass-tMass_);
423  if ( dmNew < dm )
424  {
425  dm = dmNew;
426  bestB1Idx = b1Idx;
427  bestB2Idx = b2Idx;
428  }
429  }
430  }
431  if ( bestB1Idx < 0 ) break;
432  const auto& bJet1 = jets->at(bestB1Idx);
433  const auto& bJet2 = jets->at(bestB2Idx);
434  const LorentzVector t1LVec = w1LVec + bJet1.p4();
435  const LorentzVector t2LVec = w2LVec + bJet2.p4();
436 
437  // Put all of them into candidate collection
438  if ( true ) // Trick to restrict variables' scope to avoid collision
439  {
440  const int lep1Q = lepton1.charge();
441  const int lep2Q = lepton2.charge();
442 
443  reco::GenParticle t1(lep1Q*2/3, t1LVec, genVertex_, lep1Q*6, 3, true);
444  reco::GenParticle w1(lep1Q, w1LVec, genVertex_, lep1Q*24, 3, true);
445  reco::GenParticle b1(0, bJet1.p4(), genVertex_, lep1Q*5, 1, true);
446  reco::GenParticle l1(lep1Q, lepton1.p4(), genVertex_, lepton1.pdgId(), 1, true);
447  reco::GenParticle n1(0, neutrino1.p4(), genVertex_, neutrino1.pdgId(), 1, true);
448 
449  reco::GenParticle t2(lep2Q*2/3, t2LVec, genVertex_, lep2Q*6, 3, true);
450  reco::GenParticle w2(lep2Q, w2LVec, genVertex_, lep2Q*24, 3, true);
451  reco::GenParticle b2(0, bJet2.p4(), genVertex_, lep2Q*5, 1, true);
452  reco::GenParticle l2(0, lepton2.p4(), genVertex_, lepton2.pdgId(), 1, true);
453  reco::GenParticle n2(0, neutrino2.p4(), genVertex_, neutrino2.pdgId(), 1, true);
454 
455  pseudoTop->push_back(t1);
456  pseudoTop->push_back(w1);
457  pseudoTop->push_back(b1);
458  pseudoTop->push_back(l1);
459  pseudoTop->push_back(n1);
460 
461  pseudoTop->push_back(t2);
462  pseudoTop->push_back(w2);
463  pseudoTop->push_back(b2);
464  pseudoTop->push_back(l2);
465  pseudoTop->push_back(n2);
466  }
467  }
468 
469  if ( pseudoTop->size() == 10 ) // If pseudtop decay tree is completed
470  {
471  // t->W+b
472  pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 1)); // t->W
473  pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 2)); // t->b
474  pseudoTop->at(1).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->W
475  pseudoTop->at(2).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->b
476 
477  // W->lv or W->jj
478  pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 3));
479  pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 4));
480  pseudoTop->at(3).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
481  pseudoTop->at(4).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
482 
483  // tbar->W-b
484  pseudoTop->at(5).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 6));
485  pseudoTop->at(5).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 7));
486  pseudoTop->at(6).addMother(reco::GenParticleRef(pseudoTopRefHandle, 5));
487  pseudoTop->at(7).addMother(reco::GenParticleRef(pseudoTopRefHandle, 5));
488 
489  // W->jj
490  pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 8));
491  pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 9));
492  pseudoTop->at(8).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
493  pseudoTop->at(9).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
494  }
495  } while (false);
496 
497  event.put(neutrinos, "neutrinos");
498  event.put(leptons, "leptons");
499  event.put(jets, "jets");
500 
501  event.put(pseudoTop);
502 }
bool isBHadron(const reco::Candidate *p) const
const double leptonMaxEta_
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
int i
Definition: DBlmapReader.cc:9
reco::Particle::Point genVertex_
virtual double energy() const =0
energy
const double leptonMinPt_
const double jetMaxEta_
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
edm::EDGetTokenT< edm::View< reco::Candidate > > finalStateToken_
tuple bJet1
do the razor with one or two b jets (medium CSV)
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
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::EventIDconst &, edm::Timestampconst & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
virtual double pt() const =0
transverse momentum
std::shared_ptr< JetDef > fjJetDef_
virtual int status() const =0
status word
std::vector< GenJet > GenJetCollection
collection of GenJet objects
virtual double pz() const =0
z coordinate of momentum vector
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
bool isFromHadron(const reco::Candidate *p) const
math::XYZTLorentzVector LorentzVector
std::shared_ptr< JetDef > fjLepDef_
virtual double p() const =0
magnitude of momentum vector
auto const T2 &decltype(t1.eta()) t2
Definition: deltaR.h:16
bool isnan(float x)
Definition: math.h:13
vector< PseudoJet > jets
tuple dm
Definition: symbols.py:65
bool isNull() const
Checks for null.
Definition: Ptr.h:165
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
Jets made from MC generator particles.
Definition: GenJet.h:24
virtual const Point & vertex() const =0
vertex position
virtual int charge() const =0
electric charge
virtual double py() const =0
y coordinate of momentum vector
edm::EDGetTokenT< edm::View< reco::Candidate > > genParticleToken_
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:169
virtual double px() const =0
x coordinate of momentum vector
virtual int pdgId() const =0
PDG identifier.
const double jetMinPt_
tuple status
Definition: mps_update.py:57
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
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector

Member Data Documentation

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

Definition at line 31 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

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

Definition at line 37 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

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

Definition at line 37 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

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

Definition at line 32 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

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

Definition at line 38 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

const double PseudoTopProducer::jetMaxEta_
private

Definition at line 33 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::jetMinPt_
private

Definition at line 33 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::leptonMaxEta_
private

Definition at line 33 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::leptonMinPt_
private

Definition at line 33 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::tMass_
private

Definition at line 34 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::wMass_
private

Definition at line 34 of file PseudoTopProducer.h.

Referenced by produce().