#include <PseudoTopProducer.h>

Inheritance diagram for PseudoTopProducer:

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

Public Member Functions inherited from edm::stream::EDProducerBase
	EDProducerBase ()

ModuleDescription const &	moduleDescription () const

virtual	~EDProducerBase ()

Public Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

std::vector < edm::ProductResolverIndex > const &	indiciesForPutProducts (BranchType iBranchType) const

	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

void	resolvePutIndicies (BranchType iBranchType, std::unordered_multimap< std::string, edm::ProductResolverIndex > const &iIndicies, std::string const &moduleLabel)

virtual	~ProducerBase () noexcept(false)

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector < ProductResolverIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

virtual	~EDConsumerBase () noexcept(false)

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

Public Types inherited from edm::stream::EDProducerBase
typedef EDProducerAdaptorBase	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::stream::EDProducerBase
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 15 of file PseudoTopProducer.h.

Member Typedef Documentation

typedef fastjet::JetDefinition PseudoTopProducer::JetDef

private

Definition at line 43 of file PseudoTopProducer.h.

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

private

Definition at line 30 of file PseudoTopProducer.h.

Constructor & Destructor Documentation

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

Definition at line 13 of file PseudoTopProducer.cc.

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

                                                                :
   finalStateToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("finalStates"))),
   genParticleToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("genParticles"))),
   minLeptonPt_(pset.getParameter<double>("minLeptonPt")),
   maxLeptonEta_(pset.getParameter<double>("maxLeptonEta")),
   minJetPt_(pset.getParameter<double>("minJetPt")),
   maxJetEta_(pset.getParameter<double>("maxJetEta")),
   wMass_(pset.getParameter<double>("wMass")),
   tMass_(pset.getParameter<double>("tMass")),
   minLeptonPtDilepton_(pset.getParameter<double>("minLeptonPtDilepton")),
   maxLeptonEtaDilepton_(pset.getParameter<double>("maxLeptonEtaDilepton")),
   minDileptonMassDilepton_(pset.getParameter<double>("minDileptonMassDilepton")),
   minLeptonPtSemilepton_(pset.getParameter<double>("minLeptonPtSemilepton")),
   maxLeptonEtaSemilepton_(pset.getParameter<double>("maxLeptonEtaSemilepton")),
   minVetoLeptonPtSemilepton_(pset.getParameter<double>("minVetoLeptonPtSemilepton")),
   maxVetoLeptonEtaSemilepton_(pset.getParameter<double>("maxVetoLeptonEtaSemilepton")),
   minMETSemiLepton_(pset.getParameter<double>("minMETSemiLepton")),
   minMtWSemiLepton_(pset.getParameter<double>("minMtWSemiLepton"))
 {
   const double leptonConeSize = pset.getParameter<double>("leptonConeSize");
   const double jetConeSize = pset.getParameter<double>("jetConeSize");
   fjLepDef_ = std::shared_ptr<JetDef>(new JetDef(fastjet::antikt_algorithm, leptonConeSize));
   fjJetDef_ = std::shared_ptr<JetDef>(new JetDef(fastjet::antikt_algorithm, jetConeSize));
 
   genVertex_ = reco::Particle::Point(0,0,0);
 
   produces<reco::GenParticleCollection>("neutrinos");
   produces<reco::GenJetCollection>("leptons");
   produces<reco::GenJetCollection>("jets");
 
   produces<reco::GenParticleCollection>();
 
 }

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

 {
   reco::GenParticle pOut(*dynamic_cast<const reco::GenParticle*>(p));
   pOut.clearMothers();
   pOut.clearDaughters();
   pOut.resetMothers(refHandle.id());
   pOut.resetDaughters(refHandle.id());
 
   outColl->push_back(pOut);
 
   return reco::GenParticleRef(refHandle, outColl->size()-1);
 }

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

private

Definition at line 448 of file PseudoTopProducer.cc.

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

 {
   for ( size_t i=0, n=p->numberOfDaughters(); i<n; ++i )
   {
     const reco::Candidate* dau = p->daughter(i);
     if ( p->pdgId() == dau->pdgId() ) return getLast(dau);
   }
   return p;
 }

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

private

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

private

Definition at line 473 of file PseudoTopProducer.cc.

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

Referenced by produce().

 {
   const unsigned int absPdgId = abs(p->pdgId());
   if ( !isBHadron(absPdgId) ) return false;
 
   // Do not consider this particle if it has B hadron daughter
   // For example, B* -> B0 + photon; then we drop B* and take B0 only
   for ( int i=0, n=p->numberOfDaughters(); i<n; ++i )
   {
     const reco::Candidate* dau = p->daughter(i);
     if ( isBHadron(abs(dau->pdgId())) ) return false;
   }
 
   return true;
 }

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

private

Definition at line 489 of file PseudoTopProducer.cc.

 {
   if ( absPdgId <= 100 ) return false; // Fundamental particles and MC internals
   if ( absPdgId >= 1000000000 ) return false; // Nuclei, +-10LZZZAAAI
 
   // General form of PDG ID is 7 digit form
   // +- n nr nL nq1 nq2 nq3 nJ
   //const int nJ = absPdgId % 10; // Spin
   const int nq3 = (absPdgId / 10) % 10;
   const int nq2 = (absPdgId / 100) % 10;
   const int nq1 = (absPdgId / 1000) % 10;
 
   if ( nq3 == 0 ) return false; // Diquarks
   if ( nq1 == 0 and nq2 == 5 ) return true; // B mesons
   if ( nq1 == 5 ) return true; // B baryons
 
   return false;
 }

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

private

Definition at line 458 of file PseudoTopProducer.cc.

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

Referenced by produce().

 {
   for ( size_t i=0, n=p->numberOfMothers(); i<n; ++i )
   {
     const reco::Candidate* mother = p->mother(i);
     if ( !mother ) continue;
     if ( mother->status() == 4 or mother->numberOfMothers() == 0 ) continue; // Skip incident beam
     const int pdgId = abs(mother->pdgId());
 
     if ( pdgId > 100 ) return true;
     else if ( isFromHadron(mother) ) return true;
   }
   return false;
 }

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

overridevirtual

Build jets

Implements edm::stream::EDProducerBase.

Definition at line 47 of file PseudoTopProducer.cc.

References funct::abs(), reco::CompositeRefCandidateT< D >::addDaughter(), RazorAnalyzer::bJet1, RazorAnalyzer::bJet2, reco::Candidate::charge(), funct::cos(), deltaR(), symbols::dm, PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, 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, fwrapper::jets, HLT_25ns10e33_v2_cff::leptons, visualization-live-secondInstance_cfg::m, ResonanceBuilder::mass, maxJetEta_, maxLeptonEta_, maxLeptonEtaDilepton_, maxLeptonEtaSemilepton_, maxVetoLeptonEtaSemilepton_, minDileptonMassDilepton_, HLT_25ns10e33_v2_cff::minDR, minJetPt_, minLeptonPt_, minLeptonPtDilepton_, minLeptonPtSemilepton_, minMETSemiLepton_, minMtWSemiLepton_, minVetoLeptonPtSemilepton_, reco::Candidate::mother(), eostools::move(), gen::n, GetRecoTauVFromDQM_MC_cff::next, 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::LeafCandidate::px(), reco::Candidate::py(), reco::Candidate::pz(), lumiQueryAPI::q, q1, q2, pileupReCalc_HLTpaths::scale, mathSSE::sqrt(), mps_update::status, reco::Candidate::status(), reco::t2, tMass_, MetAnalyzer::u2, reco::Candidate::vertex(), w2, wMass_, and reco::writeSpecific().

 {
   edm::Handle<edm::View<reco::Candidate> > finalStateHandle;
   event.getByToken(finalStateToken_, finalStateHandle);
 
   edm::Handle<edm::View<reco::Candidate> > genParticleHandle;
   event.getByToken(genParticleToken_, genParticleHandle);
 
   std::unique_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
   std::unique_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
   std::unique_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
   auto neutrinosRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("neutrinos");
   auto leptonsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("leptons");
   auto jetsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("jets");
 
   std::unique_ptr<reco::GenParticleCollection> pseudoTop(new reco::GenParticleCollection);
   auto pseudoTopRefHandle = event.getRefBeforePut<reco::GenParticleCollection>();
 
   // Collect unstable B-hadrons
   std::set<int> bHadronIdxs;
   for ( size_t i=0, n=genParticleHandle->size(); i<n; ++i ) {
     const reco::Candidate& p = genParticleHandle->at(i);
     const int status = p.status();
     if ( status == 1 ) continue;
 
     // Collect B-hadrons, to be used in b tagging
     if ( isBHadron(&p) ) bHadronIdxs.insert(-i-1);
   }
 
   // Collect stable leptons and neutrinos
   size_t nStables = 0;
   std::vector<size_t> leptonIdxs;
   for ( size_t i=0, n=finalStateHandle->size(); i<n; ++i ) {
     const reco::Candidate& p = finalStateHandle->at(i);
     const int absPdgId = abs(p.pdgId());
     if ( p.status() != 1 ) continue;
 
     ++nStables;
     if ( p.numberOfMothers() == 0 ) continue; // Skip orphans (if exists)
     if ( p.mother()->status() == 4 ) continue; // Treat particle as hadronic if directly from the incident beam (protect orphans in MINIAOD)
     if ( isFromHadron(&p) ) continue;
     switch ( absPdgId ) {
       case 11: case 13: // Leptons
       case 22: // Photons
         leptonIdxs.push_back(i);
         break;
       case 12: case 14: case 16:
         neutrinos->push_back(reco::GenParticle(p.charge(), p.p4(), p.vertex(), p.pdgId(), p.status(), true));
         break;
     }
   }
 
   // Sort neutrinos by pT.
   std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
 
   // Make dressed leptons with anti-kt(0.1) algorithm
   std::vector<fastjet::PseudoJet> fjLepInputs;
   fjLepInputs.reserve(leptonIdxs.size());
   for ( auto index : leptonIdxs ) {
     const reco::Candidate& p = finalStateHandle->at(index);
     if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
 
     fjLepInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
     fjLepInputs.back().set_user_index(index);
   }
 
   fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *fjLepDef_);
   std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(minLeptonPt_));
 
   leptons->reserve(fjLepJets.size());
   std::set<size_t> lepDauIdxs; // keep lepton constituents to remove from GenJet construction
   for ( auto& fjJet : fjLepJets ) {
     if ( abs(fjJet.eta()) > maxLeptonEta_ ) continue;
 
     // Get jet constituents from fastJet
     const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
     // Convert to CandidatePtr
     std::vector<reco::CandidatePtr> constituents;
     reco::CandidatePtr lepCand;
     for ( auto& fjConstituent : fjConstituents ) {
       const size_t index = fjConstituent.user_index();
       reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
       const int absPdgId = abs(cand->pdgId());
       if ( absPdgId == 11 or absPdgId == 13 ) {
         if ( lepCand.isNonnull() and lepCand->pt() > cand->pt() ) continue; // Choose one with highest pt
         lepCand = cand;
       }
       constituents.push_back(cand);
     }
     if ( lepCand.isNull() ) continue;
 
     const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
     reco::GenJet lepJet;
     reco::writeSpecific(lepJet, jetP4, genVertex_, constituents, eventSetup);
 
     lepJet.setPdgId(lepCand->pdgId());
     lepJet.setCharge(lepCand->charge());
 
     const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
     lepJet.setJetArea(jetArea);
 
     leptons->push_back(lepJet);
 
     // Keep constituent indices to be used in the next step.
     for ( auto& fjConstituent : fjConstituents ) {
       lepDauIdxs.insert(fjConstituent.user_index());
     }
   }
 
   // Now proceed to jets.
   // Jets: anti-kt excluding the e, mu, nu, and photons in selected leptons.
   std::vector<fastjet::PseudoJet> fjJetInputs;
   fjJetInputs.reserve(nStables);
   for ( size_t i=0, n=finalStateHandle->size(); i<n; ++i ) {
     const reco::Candidate& p = finalStateHandle->at(i);
     if ( p.status() != 1 ) continue;
     if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
 
     const int absId = std::abs(p.pdgId());
     if ( absId == 12 or absId == 14 or absId == 16 ) continue;
     if ( lepDauIdxs.find(i) != lepDauIdxs.end() ) continue;
 
     fjJetInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
     fjJetInputs.back().set_user_index(i);
   }
   for ( auto index : bHadronIdxs ) {
     const reco::Candidate& p = genParticleHandle->at(abs(index+1));
     if ( std::isnan(p.pt()) or p.pt() <= 0 ) continue;
 
     const double scale = 1e-20/p.p();
     fjJetInputs.push_back(fastjet::PseudoJet(p.px()*scale, p.py()*scale, p.pz()*scale, p.energy()*scale));
     fjJetInputs.back().set_user_index(index);
   }
 
   fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *fjJetDef_);
   std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(minJetPt_));
 
   jets->reserve(fjJets.size());
   std::vector<size_t> bjetIdxs, ljetIdxs;
   for ( auto& fjJet : fjJets ) {
     if ( abs(fjJet.eta()) > maxJetEta_ ) continue;
 
     // Get jet constituents from fastJet
     const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
     // Convert to CandidatePtr
     std::vector<reco::CandidatePtr> constituents;
     bool hasBHadron = false;
     for ( size_t j=0, m=fjConstituents.size(); j<m; ++j ) {
       const int index = fjConstituents[j].user_index();
       if ( bHadronIdxs.find(index) != bHadronIdxs.end() ) {
         hasBHadron = true;
         continue;
       }
       reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
       constituents.push_back(cand);
     }
 
     const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
     reco::GenJet genJet;
     reco::writeSpecific(genJet, jetP4, genVertex_, constituents, eventSetup);
 
     const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
     genJet.setJetArea(jetArea);
     if ( hasBHadron ) {
       genJet.setPdgId(5);
       bjetIdxs.push_back(jets->size());
     }
     else {
       ljetIdxs.push_back(jets->size());
     }
 
     jets->push_back(genJet);
   }
 
   // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
   // NOTE : A C++ trick, use do-while instead of long-nested if-statements.
   do {
     if ( bjetIdxs.size() < 2 ) break;
 
     // Note : we will do dilepton or semilepton channel only
     if ( leptons->size() == 2 and neutrinos->size() >= 2 ) {
       // Start from dilepton channel
       const int q1 = leptons->at(0).charge();
       const int q2 = leptons->at(1).charge();
       if ( q1*q2 > 0 ) break;
 
       const auto& lepton1 = q1 > 0 ? leptons->at(0) : leptons->at(1);
       const auto& lepton2 = q1 > 0 ? leptons->at(1) : leptons->at(0);
 
       if ( lepton1.pt() < minLeptonPtDilepton_ or std::abs(lepton1.eta()) > maxLeptonEtaDilepton_ ) break;
       if ( lepton2.pt() < minLeptonPtDilepton_ or std::abs(lepton2.eta()) > maxLeptonEtaDilepton_ ) break;
       if ( (lepton1.p4()+lepton2.p4()).mass() < minDileptonMassDilepton_ ) break;
 
       double dm = 1e9;
       int selNu1 = -1, selNu2 = -1;
       for ( int i=0; i<2; ++i ) { // Consider leading 2 neutrinos. Neutrino vector is already sorted by pT
         const double dm1 = std::abs((lepton1.p4()+neutrinos->at(i).p4()).mass()-wMass_);
         for ( int j=0; j<2; ++j ) {
           if ( i == j ) continue;
           const double dm2 = std::abs((lepton2.p4()+neutrinos->at(j).p4()).mass()-wMass_);
           const double newDm = dm1+dm2;
 
           if ( newDm < dm ) { dm = newDm; selNu1 = i; selNu2 = j; }
         }
       }
       if ( dm >= 1e9 ) break;
 
       const auto& nu1 = neutrinos->at(selNu1);
       const auto& nu2 = neutrinos->at(selNu2);
       const auto w1LVec = lepton1.p4()+nu1.p4();
       const auto w2LVec = lepton2.p4()+nu2.p4();
 
       // Contiue to top quarks
       dm = 1e9; // Reset once again for top combination.
       int selB1 = -1, selB2 = -1;
       for ( auto i : bjetIdxs ) {
         const double dm1 = std::abs((w1LVec+jets->at(i).p4()).mass()-tMass_);
         for ( auto j : bjetIdxs ) {
           if ( i == j ) continue;
           const double dm2 = std::abs((w2LVec+jets->at(j).p4()).mass()-tMass_);
           const double newDm = dm1+dm2;
 
           if ( newDm < dm ) { dm = newDm; selB1 = i; selB2 = j; }
         }
       }
       if ( dm >= 1e9 ) break;
 
       const auto& bJet1 = jets->at(selB1);
       const auto& bJet2 = jets->at(selB2);
       const auto t1LVec = w1LVec + bJet1.p4();
       const auto t2LVec = w2LVec + bJet2.p4();
 
       // Recalculate lepton charges (needed due to initialization of leptons wrt sign of a charge)
       const int lepQ1 = lepton1.charge();
       const int lepQ2 = lepton2.charge();
 
       // Put all of them into candidate collection
       reco::GenParticle t1(lepQ1*2/3., t1LVec, genVertex_, lepQ1*6, 3, false);
       reco::GenParticle w1(lepQ1, w1LVec, genVertex_, lepQ1*24, 3, true);
       reco::GenParticle b1(-lepQ1/3., bJet1.p4(), genVertex_, lepQ1*5, 1, true);
       reco::GenParticle l1(lepQ1, lepton1.p4(), genVertex_, lepton1.pdgId(), 1, true);
       reco::GenParticle n1(0, nu1.p4(), genVertex_, nu1.pdgId(), 1, true);
 
       reco::GenParticle t2(lepQ2*2/3., t2LVec, genVertex_, lepQ2*6, 3, false);
       reco::GenParticle w2(lepQ2, w2LVec, genVertex_, lepQ2*24, 3, true);
       reco::GenParticle b2(-lepQ2/3., bJet2.p4(), genVertex_, lepQ2*5, 1, true);
       reco::GenParticle l2(lepQ2, lepton2.p4(), genVertex_, lepton2.pdgId(), 1, true);
       reco::GenParticle n2(0, nu2.p4(), genVertex_, nu2.pdgId(), 1, true);
 
       pseudoTop->push_back(t1);
       pseudoTop->push_back(t2);
 
       pseudoTop->push_back(w1);
       pseudoTop->push_back(b1);
       pseudoTop->push_back(l1);
       pseudoTop->push_back(n1);
 
       pseudoTop->push_back(w2);
       pseudoTop->push_back(b2);
       pseudoTop->push_back(l2);
       pseudoTop->push_back(n2);
     }
     else if ( leptons->size() >= 1 and neutrinos->size() >= 1 and ljetIdxs.size() >= 2 ) {
       // Then continue to the semi-leptonic channel
       const auto& lepton = leptons->at(0);
       if ( lepton.pt() < minLeptonPtSemilepton_ or std::abs(lepton.eta()) > maxLeptonEtaSemilepton_ ) break;
 
       // Skip event if there are veto leptons
       bool hasVetoLepton = false;
       for ( auto vetoLeptonItr = next(leptons->begin()); vetoLeptonItr != leptons->end(); ++vetoLeptonItr ) {
         if ( vetoLeptonItr->pt() > minVetoLeptonPtSemilepton_ and
              std::abs(vetoLeptonItr->eta()) < maxVetoLeptonEtaSemilepton_ ) {
           hasVetoLepton = true;
         }
       }
       if ( hasVetoLepton ) break;
 
       // Calculate MET
       double metX = 0, metY = 0;
       for ( auto neutrino : *neutrinos ) {
         metX += neutrino.px();
         metY += neutrino.py();
       }
       const double metPt = std::hypot(metX, metY);
       if ( metPt < minMETSemiLepton_ ) break;
 
       const double mtW = std::sqrt( 2*lepton.pt()*metPt*cos(lepton.phi()-atan2(metX, metY)) );
       if ( mtW < minMtWSemiLepton_ ) break;
 
       // Solve pz to give wMass_^2 = (lepE+energy)^2 - (lepPx+metX)^2 - (lepPy+metY)^2 - (lepPz+pz)^2
       // -> (wMass_^2)/2 = lepE*sqrt(metPt^2+pz^2) - lepPx*metX - lepPy*metY - lepPz*pz
       // -> C := (wMass_^2)/2 + lepPx*metX + lepPy*metY
       // -> lepE^2*(metPt^2+pz^2) = C^2 + 2*C*lepPz*pz + lepPz^2*pz^2
       // -> (lepE^2-lepPz^2)*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
       // -> lepPt^2*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
       const double lpz = lepton.pz(), le = lepton.energy(), lpt = lepton.pt();
       const double cf = (wMass_*wMass_)/2 + lepton.px()*metX + lepton.py()*metY;
       const double det = cf*cf*lpz*lpz - lpt*lpt*(le*le*metPt*metPt - cf*cf);
       const double pz = (cf*lpz + (cf < 0 ? -sqrt(det) : sqrt(det)))/lpt/lpt;
       const reco::Candidate::LorentzVector nu1P4(metX, metY, pz, std::hypot(metPt, pz));
       const auto w1LVec = lepton.p4()+nu1P4;
 
       // Continue to build leptonic pseudo top
       double minDR = 1e9;
       int selB1 = -1;
       for ( auto b1Itr = bjetIdxs.begin(); b1Itr != bjetIdxs.end(); ++b1Itr ) {
         const double dR = deltaR(jets->at(*b1Itr).p4(), w1LVec);
         if ( dR < minDR ) { selB1 = *b1Itr; minDR = dR; }
       }
       if ( selB1 == -1 ) break;
       const auto& bJet1 = jets->at(selB1);
       const auto t1LVec = w1LVec + bJet1.p4();
 
       // Build hadronic pseudo W, take leading two jets (ljetIdxs are (implicitly) sorted by pT)
       const auto& wJet1 = jets->at(ljetIdxs[0]);
       const auto& wJet2 = jets->at(ljetIdxs[1]);
       const auto w2LVec = wJet1.p4() + wJet2.p4();
 
       // Contiue to top quarks
       double dm = 1e9;
       int selB2 = -1;
       for ( auto i : bjetIdxs ) {
         if ( int(i) == selB1 ) continue;
         const double newDm = std::abs((w2LVec+jets->at(i).p4()).mass()-tMass_);
         if ( newDm < dm ) { dm = newDm; selB2 = i; }
       }
       if ( dm >= 1e9 ) break;
 
       const auto& bJet2 = jets->at(selB2);
       const auto t2LVec = w2LVec + bJet2.p4();
 
       const int q = lepton.charge();
       // Put all of them into candidate collection
       reco::GenParticle t1(q*2/3., t1LVec, genVertex_, q*6, 3, false);
       reco::GenParticle w1(q, w1LVec, genVertex_, q*24, 3, true);
       reco::GenParticle b1(-q/3., bJet1.p4(), genVertex_, q*5, 1, true);
       reco::GenParticle l1(q, lepton.p4(), genVertex_, lepton.pdgId(), 1, true);
       reco::GenParticle n1(0, nu1P4, genVertex_, q*(std::abs(lepton.pdgId())+1), 1, true);
 
       reco::GenParticle t2(-q*2/3., t2LVec, genVertex_, -q*6, 3, false);
       reco::GenParticle w2(-q, w2LVec, genVertex_, -q*24, 3, true);
       reco::GenParticle b2(0, bJet2.p4(), genVertex_, -q*5, 1, true);
       reco::GenParticle u2(0, wJet1.p4(), genVertex_, -2*q, 1, true);
       reco::GenParticle d2(0, wJet2.p4(), genVertex_, q, 1, true);
 
       pseudoTop->push_back(t1);
       pseudoTop->push_back(t2);
 
       pseudoTop->push_back(w1);
       pseudoTop->push_back(b1);
       pseudoTop->push_back(l1);
       pseudoTop->push_back(n1);
 
       pseudoTop->push_back(w2);
       pseudoTop->push_back(b2);
       pseudoTop->push_back(u2);
       pseudoTop->push_back(d2);
     }
   } while (false);
 
   if ( pseudoTop->size() == 10 ) // If pseudtop decay tree is completed
   {
     // t->W+b
     pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 2)); // t->W
     pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 3)); // t->b
     pseudoTop->at(2).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->W
     pseudoTop->at(3).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0)); // t->b
 
     // W->lv or W->jj
     pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 4));
     pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 5));
     pseudoTop->at(4).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
     pseudoTop->at(5).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
 
     // tbar->W-b
     pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 6));
     pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 7));
     pseudoTop->at(6).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
     pseudoTop->at(7).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
 
     // W->jj
     pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 8));
     pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 9));
     pseudoTop->at(8).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
     pseudoTop->at(9).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
   }
 
   event.put(std::move(neutrinos), "neutrinos");
   event.put(std::move(leptons), "leptons");
   event.put(std::move(jets), "jets");
 
   event.put(std::move(pseudoTop));
 }

Member Data Documentation

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

private

Definition at line 33 of file PseudoTopProducer.h.

Referenced by produce().

std::shared_ptr<JetDef> PseudoTopProducer::fjJetDef_

private

Definition at line 44 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

std::shared_ptr<JetDef> PseudoTopProducer::fjLepDef_

private

Definition at line 44 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

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

private

Definition at line 34 of file PseudoTopProducer.h.

Referenced by produce().

reco::Particle::Point PseudoTopProducer::genVertex_

private

Definition at line 45 of file PseudoTopProducer.h.

Referenced by produce(), and PseudoTopProducer().

const double PseudoTopProducer::maxJetEta_

private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::maxLeptonEta_

private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::maxLeptonEtaDilepton_

private

Definition at line 37 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::maxLeptonEtaSemilepton_

private

Definition at line 39 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::maxVetoLeptonEtaSemilepton_

private

Definition at line 40 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minDileptonMassDilepton_

private

Definition at line 38 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minJetPt_

private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minLeptonPt_

private

Definition at line 35 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minLeptonPtDilepton_

private

Definition at line 37 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minLeptonPtSemilepton_

private

Definition at line 39 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minMETSemiLepton_

private

Definition at line 41 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minMtWSemiLepton_

private

Definition at line 41 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::minVetoLeptonPtSemilepton_

private

Definition at line 40 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::tMass_

private

Definition at line 36 of file PseudoTopProducer.h.

Referenced by produce().

const double PseudoTopProducer::wMass_

private

Definition at line 36 of file PseudoTopProducer.h.

Referenced by produce().

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation