Module to produce the subset of generator particles directly contained in top decay chains. More...

#include "TopQuarkAnalysis/TopEventProducers/interface/TopDecaySubset.h"

Inheritance diagram for TopDecaySubset:

Public Types
enum	FillMode { kStable, kME }

enum	RunMode { kRun1, kRun2 }
	supported modes to run the code More...

enum	ShowerModel { kStart = -1, kNone, kPythia, kHerwig, kPythia8, kSherpa }
	classification of potential shower types More...

Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

Public Types inherited from edm::ProducerBase
using	ModuleToResolverIndicies = std::unordered_multimap< std::string, std::tuple< edm::TypeID const , const char , edm::ProductResolverIndex > >

typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Public Member Functions
void	produce (edm::Event &event, const edm::EventSetup &setup) override
	write output into the event More...

	TopDecaySubset (const edm::ParameterSet &cfg)
	default constructor More...

	~TopDecaySubset () override
	default destructor More...

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

SerialTaskQueue *	globalLuminosityBlocksQueue ()

SerialTaskQueue *	globalRunsQueue ()

ModuleDescription const &	moduleDescription () const

	~EDProducer () override

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

std::vector< edm::ProductResolverIndex > const &	putTokenIndexToProductResolverIndex () const

std::vector< bool > const &	recordProvenanceList () const

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, ModuleToResolverIndicies const &iIndicies, std::string const &moduleLabel)

TypeLabelList const &	typeLabelList () const
	used by the fwk to register the list of products of this module More...

	~ProducerBase () noexcept(false) override

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

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase &&)=default

	EDConsumerBase (EDConsumerBase const &)=delete

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

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 &	itemsToGetFrom (BranchType iType) 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 &	operator= (EDConsumerBase &&)=default

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

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

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

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
void	addDaughters (int &idx, const reco::GenParticle::const_iterator part, reco::GenParticleCollection &target, bool recursive=true)
	recursively fill vector for all further decay particles of a given particle More...

void	addRadiation (int &idx, const reco::GenParticle *part, reco::GenParticleCollection &target)

void	addRadiation (int &idx, const reco::GenParticle::const_iterator part, reco::GenParticleCollection &target)
	fill vector including all radiations from quarks originating from W/top More...

ShowerModel	checkShowerModel (const std::vector< const reco::GenParticle * > &tops) const
	check the decay chain for the used shower model More...

ShowerModel	checkShowerModel (edm::Event &event)
	check the embedded MC information for the shower model More...

void	checkWBosons (std::vector< const reco::GenParticle * > &tops) const
	check whether W bosons are contained in the original gen particle listing More...

void	clearReferences ()
	clear references More...

void	fillListing (const std::vector< const reco::GenParticle * > &primalTops, const std::vector< const reco::GenParticle * > &decayingTops, reco::GenParticleCollection &target)
	fill output vector for full decay chain More...

void	fillListing (const std::vector< const reco::GenParticle * > &tops, reco::GenParticleCollection &target)
	fill output vector for full decay chain More...

void	fillReferences (const reco::GenParticleRefProd &refProd, reco::GenParticleCollection &target)
	fill references for output vector More...

std::vector< const reco::GenParticle * >	findDecayingTops (const reco::GenParticleCollection &parts)

const reco::GenParticle *	findLastParticleInChain (const reco::GenParticle *p)

std::vector< const reco::GenParticle * >	findPrimalTops (const reco::GenParticleCollection &parts)

const reco::GenParticle *	findPrimalW (const reco::GenParticle *top)

std::vector< const reco::GenParticle * >	findTops (const reco::GenParticleCollection &parts)
	find top quarks in list of input particles More...

reco::Particle::LorentzVector	p4 (const reco::GenParticle::const_iterator part, int statusFlag)
	calculate lorentz vector from input More...

reco::Particle::LorentzVector	p4 (const std::vector< const reco::GenParticle * >::const_iterator top, int statusFlag)
	calculate lorentz vector from input (dedicated to top reconstruction) More...

Private Attributes
bool	addRadiation_
	add radiation or not? More...

FillMode	fillMode_

edm::EDGetTokenT< GenEventInfoProduct >	genEventInfo_srcToken_
	input tag for the genEventInfo source More...

int	motherPartIdx_

std::map< int, std::vector< int > >	refs_
	management of daughter indices for fillRefs More...

RunMode	runMode_
	run mode (Run1 \|\| Run2) More...

ShowerModel	showerModel_
	parton shower mode (filled in checkShowerModel) More...

edm::EDGetTokenT< reco::GenParticleCollection >	srcToken_
	input tag for the genParticle source More...

Additional Inherited Members
Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

static bool	wantsGlobalLuminosityBlocks ()

static bool	wantsGlobalRuns ()

static bool	wantsStreamLuminosityBlocks ()

static bool	wantsStreamRuns ()

Protected Member Functions inherited from edm::ProducerBase
template<class ProductType >
BranchAliasSetterT< ProductType >	produces ()
	declare what type of product will make and with which optional label More...

template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType >	produces ()

template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType >	produces ()

BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<BranchType B>
BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<Transition B>
BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

template<class ProductType >
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

ProducesCollector	producesCollector ()

Protected Member Functions inherited from edm::EDConsumerBase
EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (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 ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

Detailed Description

Module to produce the subset of generator particles directly contained in top decay chains.

The module produces the subset of generator particles directly contained in top decay chains. The particles are saved as a collection of reco::GenParticles. Depending on the configuration of the module, the 4-vector kinematics can be taken from the status-3 particles (ME before parton showering) or from the status-2 particles (after parton showering), additionally radiated gluons may be considered during the creation of the subset or not.

Definition at line 28 of file TopDecaySubset.h.

Member Enumeration Documentation

◆ FillMode

enum TopDecaySubset::FillMode

supported modes to fill the new vectors of gen particles

Enumerator
kStable
kME

Definition at line 32 of file TopDecaySubset.h.

32 { kStable, kME };

◆ RunMode

enum TopDecaySubset::RunMode

supported modes to run the code

Enumerator
kRun1
kRun2

Definition at line 36 of file TopDecaySubset.h.

36 { kRun1, kRun2 };

◆ ShowerModel

enum TopDecaySubset::ShowerModel

classification of potential shower types

Enumerator
kStart
kNone
kPythia
kHerwig
kPythia8
kSherpa

Definition at line 34 of file TopDecaySubset.h.

34 { kStart = -1, kNone, kPythia, kHerwig, kPythia8, kSherpa };

Constructor & Destructor Documentation

◆ TopDecaySubset()

TopDecaySubset::TopDecaySubset ( const edm::ParameterSet & cfg )

explicit

default constructor

Definition at line 10 of file TopDecaySubset.cc.

     : srcToken_(consumes<reco::GenParticleCollection>(cfg.getParameter<edm::InputTag>("src"))),
       genEventInfo_srcToken_(mayConsume<GenEventInfoProduct>(edm::InputTag("generator"))),
       addRadiation_(cfg.getParameter<bool>("addRadiation")),
       showerModel_(kStart),
       runMode_(kRun1) {
   // mapping of the corresponding fillMode; see FillMode
   // enumerator of TopDecaySubset for available modes
   std::string mode = cfg.getParameter<std::string>("fillMode");
   if (mode == "kME")
     fillMode_ = kME;
   else if (mode == "kStable")
     fillMode_ = kStable;
   else
     throw cms::Exception("Configuration") << mode << " is not a supported FillMode!\n";
  
   mode = cfg.getParameter<std::string>("runMode");
   if (mode == "Run1")
     runMode_ = kRun1;
   else if (mode == "Run2")
     runMode_ = kRun2;
   else
     throw cms::Exception("Configuration") << mode << " is not a supported RunMode!\n";
  
   // produces a set of GenParticles following
   // the decay branch of top quarks to the first level of
   // stable decay products
   produces<reco::GenParticleCollection>();
 }

References looper::cfg, Exception, fillMode_, kME, kRun1, kRun2, kStable, ALCARECOPromptCalibProdSiPixelAli0T_cff::mode, runMode_, and AlCaHLTBitMon_QueryRunRegistry::string.

◆ ~TopDecaySubset()

TopDecaySubset::~TopDecaySubset ( )

override

default destructor

Definition at line 41 of file TopDecaySubset.cc.

41 {}

Member Function Documentation

◆ addDaughters()

void TopDecaySubset::addDaughters	(	int &	idx,
		const reco::GenParticle::const_iterator	part,
		reco::GenParticleCollection &	target,
		bool	recursive = `true`
	)

private

recursively fill vector for all further decay particles of a given particle

Definition at line 631 of file TopDecaySubset.cc.

                                                   {
   std::vector<int> daughters;
   int idxBuffer = idx;
   for (reco::GenParticle::const_iterator daughter = part->begin(); daughter != part->end(); ++daughter) {
     std::unique_ptr<reco::GenParticle> ptr(new reco::GenParticle(
         daughter->threeCharge(), daughter->p4(), daughter->vertex(), daughter->pdgId(), daughter->status(), false));
     target.push_back(*ptr);
     // increment & push index of daughter
     daughters.push_back(++idx);
     // continue recursively if desired
     if (recursive) {
       addDaughters(idx, daughter, target);
     }
   }
   if (!daughters.empty()) {
     refs_[idxBuffer] = daughters;
   }
 }

References BPhysicsValidation_cfi::daughters, training_settings::idx, postprocess-scan-build::recursive, refs_, and filterCSVwithJSON::target.

Referenced by fillListing().

◆ addRadiation() [1/2]

void TopDecaySubset::addRadiation	(	int &	idx,
		const reco::GenParticle *	part,
		reco::GenParticleCollection &	target
	)

private

Definition at line 612 of file TopDecaySubset.cc.

                                                                                                           {
   std::vector<int> daughters;
   int idxBuffer = idx;
   for (reco::GenParticle::const_iterator daughter = part->begin(); daughter != part->end(); ++daughter) {
     // either pick daughters as long as they are different
     // to the initial particle
     if (daughter->pdgId() != part->pdgId()) {
       std::unique_ptr<reco::GenParticle> ptr(new reco::GenParticle(
           daughter->threeCharge(), daughter->p4(), daughter->vertex(), daughter->pdgId(), daughter->status(), false));
       target.push_back(*ptr);
       daughters.push_back(++idx);  //push index of daughter
     }
   }
   if (!daughters.empty()) {
     refs_[idxBuffer] = daughters;
   }
 }

References BPhysicsValidation_cfi::daughters, training_settings::idx, refs_, and filterCSVwithJSON::target.

◆ addRadiation() [2/2]

void TopDecaySubset::addRadiation	(	int &	idx,
		const reco::GenParticle::const_iterator	part,
		reco::GenParticleCollection &	target
	)

private

fill vector including all radiations from quarks originating from W/top

Definition at line 591 of file TopDecaySubset.cc.

                                                                      {
   std::vector<int> daughters;
   int idxBuffer = idx;
   for (reco::GenParticle::const_iterator daughter = part->begin(); daughter != part->end(); ++daughter) {
     if (daughter->status() <= TopDecayID::stable && daughter->pdgId() != part->pdgId()) {
       // skip comment lines and make sure that
       // the particle is not double counted as
       // daughter of itself
       std::unique_ptr<reco::GenParticle> ptr(new reco::GenParticle(
           daughter->threeCharge(), daughter->p4(), daughter->vertex(), daughter->pdgId(), daughter->status(), false));
       target.push_back(*ptr);
       daughters.push_back(++idx);  //push index of daughter
     }
   }
   if (!daughters.empty()) {
     refs_[idxBuffer] = daughters;
   }
 }

References BPhysicsValidation_cfi::daughters, training_settings::idx, refs_, TopDecayID::stable, and filterCSVwithJSON::target.

Referenced by fillListing().

◆ checkShowerModel() [1/2]

TopDecaySubset::ShowerModel TopDecaySubset::checkShowerModel ( const std::vector< const reco::GenParticle * > & tops ) const

private

check the decay chain for the used shower model

check the decay chain for the exploited shower model

Definition at line 202 of file TopDecaySubset.cc.

                                                                                                               {
   for (std::vector<const reco::GenParticle*>::const_iterator it = tops.begin(); it != tops.end(); ++it) {
     const reco::GenParticle* top = *it;
     // check for kHerwig type showers: here the status 3 top quark will
     // have a single status 2 top quark as daughter, which has again 3
     // or more status 2 daughters
     if (top->numberOfDaughters() == 1) {
       if (top->begin()->pdgId() == top->pdgId() && top->begin()->status() == TopDecayID::stable &&
           top->begin()->numberOfDaughters() > 1)
         return kHerwig;
     }
     // check for kPythia type showers: here the status 3 top quark will
     // have all decay products and a status 2 top quark as daughters
     // the status 2 top quark will be w/o further daughters
     if (top->numberOfDaughters() > 1) {
       bool containsWBoson = false, containsQuarkDaughter = false;
       for (reco::GenParticle::const_iterator td = top->begin(); td != top->end(); ++td) {
         if (std::abs(td->pdgId()) < TopDecayID::tID)
           containsQuarkDaughter = true;
         if (std::abs(td->pdgId()) == TopDecayID::WID)
           containsWBoson = true;
       }
       if (containsQuarkDaughter && containsWBoson)
         return kPythia;
     }
   }
   // if neither Herwig nor Pythia like
   if (tops.empty())
     edm::LogInfo("decayChain") << " Failed to find top quarks in decay chain. Will assume that this a \n"
                                << " non-top sample and produce an empty decaySubset.\n";
   else
     throw edm::Exception(edm::errors::LogicError,
                          " Can not find back any of the supported hadronization models. Models \n"
                          " which are supported are:                                            \n"
                          " Pythia  LO(+PS): Top(status 3) --> WBoson(status 3), Quark(status 3)\n"
                          " Herwig NLO(+PS): Top(status 2) --> Top(status 3) --> Top(status 2)  \n");
   return kNone;
 }

References funct::abs(), reco::Candidate::begin(), reco::Candidate::end(), Exception, kHerwig, kNone, kPythia, edm::errors::LogicError, reco::CompositeRefCandidateT< D >::numberOfDaughters(), reco::LeafCandidate::pdgId(), TopDecayID::stable, TopDecayID::tID, and TopDecayID::WID.

Referenced by produce().

◆ checkShowerModel() [2/2]

TopDecaySubset::ShowerModel TopDecaySubset::checkShowerModel ( edm::Event & event )

private

check the embedded MC information for the shower model

Definition at line 242 of file TopDecaySubset.cc.

                                                                           {
   edm::Handle<GenEventInfoProduct> genEvtInfoProduct;
   event.getByToken(genEventInfo_srcToken_, genEvtInfoProduct);
  
   std::string moduleName = "";
   if (genEvtInfoProduct.isValid()) {
     const edm::Provenance& prov = event.getProvenance(genEvtInfoProduct.id());
     moduleName = edm::moduleName(prov, event.processHistory());
   }
  
   ShowerModel shower(kStart);
   if (moduleName.find("Pythia6") != std::string::npos)
     shower = kPythia;
   else if (moduleName.find("Pythia8") != std::string::npos)
     shower = kPythia8;
   else if (moduleName.find("Herwig6") != std::string::npos)
     shower = kHerwig;
   else if (moduleName.find("ThePEG") != std::string::npos)
     // Herwig++
     shower = kHerwig;
   else if (moduleName.find("Sherpa") != std::string::npos)
     shower = kSherpa;
   else
     shower = kNone;
   return shower;
 }

References genEventInfo_srcToken_, edm::HandleBase::id(), edm::HandleBase::isValid(), kHerwig, kNone, kPythia, kPythia8, kSherpa, kStart, EcalCalibMonitorClient_cfi::moduleName, edm::moduleName(), and AlCaHLTBitMon_QueryRunRegistry::string.

◆ checkWBosons()

void TopDecaySubset::checkWBosons ( std::vector< const reco::GenParticle * > & tops ) const

private

check whether W bosons are contained in the original gen particle listing

check whether the W boson is contained in the original gen particle listing

Definition at line 269 of file TopDecaySubset.cc.

                                                                                {
   unsigned nTops = tops.size();
   for (std::vector<const reco::GenParticle*>::iterator it = tops.begin(); it != tops.end();) {
     const reco::GenParticle* top = *it;
     bool isContained = false;
     bool expectedStatus = false;
     if (showerModel_ != kPythia && top->begin() == top->end())
       throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && top->begin()==top->end()\n");
     for (reco::GenParticle::const_iterator td = ((showerModel_ == kPythia) ? top->begin() : top->begin()->begin());
          td != ((showerModel_ == kPythia) ? top->end() : top->begin()->end());
          ++td) {
       if (std::abs(td->pdgId()) == TopDecayID::WID) {
         isContained = true;
         if (((showerModel_ == kPythia) ? td->status() == TopDecayID::unfrag : td->status() == TopDecayID::stable)) {
           expectedStatus = true;
           break;
         }
       }
     }
     if (!expectedStatus) {
       it = tops.erase(it);
       if (isContained)
         edm::LogInfo("decayChain") << " W boson does not have the expected status. This happens, e.g.,      \n"
                                    << " with MC@NLO in the case of additional ttbar pairs from radiated     \n"
                                    << " gluons. We hope everything is fine, remove the correspdonding top   \n"
                                    << " quark from our list since it is not part of the primary ttbar pair  \n"
                                    << " and try to continue.                                                \n";
     } else
       it++;
   }
   if (tops.empty() && nTops != 0)
     throw edm::Exception(edm::errors::LogicError,
                          " Did not find a W boson with appropriate status for any of the top   \n"
                          " quarks in this event. This means that the hadronization of the W    \n"
                          " boson might be screwed up or there is another problem with the      \n"
                          " particle listing in this MC sample. Please contact an expert.       \n");
 }

References funct::abs(), reco::Candidate::begin(), reco::Candidate::end(), kPythia, edm::errors::LogicError, showerModel_, TopDecayID::stable, TopDecayID::unfrag, and TopDecayID::WID.

Referenced by produce().

◆ clearReferences()

void TopDecaySubset::clearReferences ( )

private

clear references

Definition at line 654 of file TopDecaySubset.cc.

                                      {
   // clear vector of references
   refs_.clear();
   // set idx for mother particles to a start value
   // of -1 (the first entry will raise it to 0)
   motherPartIdx_ = -1;
 }

References motherPartIdx_, and refs_.

Referenced by produce().

◆ fillListing() [1/2]

void TopDecaySubset::fillListing	(	const std::vector< const reco::GenParticle * > &	primalTops,
		const std::vector< const reco::GenParticle * > &	decayingTops,
		reco::GenParticleCollection &	target
	)

private

fill output vector for full decay chain

Definition at line 427 of file TopDecaySubset.cc.

                                                                     {
   std::vector<const reco::GenParticle*>::const_iterator top_start;
   std::vector<const reco::GenParticle*>::const_iterator top_end;
   if (fillMode_ == kME) {
     top_start = primalTops.begin();
     top_end = primalTops.end();
   } else {
     top_start = decayingTops.begin();
     top_end = decayingTops.end();
   }
   for (std::vector<const reco::GenParticle*>::const_iterator it = top_start; it != top_end; ++it) {
     const reco::GenParticle* t = *it;
     // summation might happen here
     std::unique_ptr<reco::GenParticle> topPtr(
         new reco::GenParticle(t->threeCharge(), t->p4(), t->vertex(), t->pdgId(), t->status(), false));
     target.push_back(*topPtr);
     ++motherPartIdx_;
  
     // keep the top index for the map to manage the daughter refs
     int iTop = motherPartIdx_;
     std::vector<int> topDaughters;
     // define the W boson index (to be set later) for the map to
     // manage the daughter refs
     int iW = 0;
     std::vector<int> wDaughters;
     const reco::GenParticle* final_top = findLastParticleInChain(t);
  
     //iterate over top daughters
     for (reco::GenParticle::const_iterator td = final_top->begin(); td != final_top->end(); ++td) {
       if (std::abs(td->pdgId()) <= TopDecayID::bID) {
         // if particle is beauty or other quark
         std::unique_ptr<reco::GenParticle> qPtr(
             new reco::GenParticle(td->threeCharge(), td->p4(), td->vertex(), td->pdgId(), td->status(), false));
         target.push_back(*qPtr);
         // increment & push index of the top daughter
         topDaughters.push_back(++motherPartIdx_);
         if (addRadiation_) {
           // for radation to be added we first need to
           // pick the last quark in the MC chain
           const reco::GenParticle* last_q = findLastParticleInChain(static_cast<const reco::GenParticle*>(&*td));
           addRadiation(motherPartIdx_, last_q, target);
         }
       } else if (std::abs(td->pdgId()) == TopDecayID::WID) {
         // ladies and gentleman, we have a W boson
         std::unique_ptr<reco::GenParticle> WPtr(
             new reco::GenParticle(td->threeCharge(), td->p4(), td->vertex(), td->pdgId(), td->status(), false));
         target.push_back(*WPtr);
         // increment & push index of the top daughter
         topDaughters.push_back(++motherPartIdx_);
         iW = motherPartIdx_;
  
         // next are the daughers of our W boson
         // for Pythia 6 this is wrong as the last W has no daughters at all!
         // instead the status 3 W has 3 daughters: q qbar' and W (WTF??!)
         const reco::GenParticle* decaying_W = findLastParticleInChain(static_cast<const reco::GenParticle*>(&*td));
         for (reco::GenParticle::const_iterator wd = decaying_W->begin(); wd != decaying_W->end(); ++wd) {
           if (!(std::abs(wd->pdgId()) == TopDecayID::WID)) {
             std::unique_ptr<reco::GenParticle> qPtr(
                 new reco::GenParticle(wd->threeCharge(), wd->p4(), wd->vertex(), wd->pdgId(), wd->status(), false));
             target.push_back(*qPtr);
             // increment & push index of the top daughter
             wDaughters.push_back(++motherPartIdx_);
             const reco::GenParticle* last_q = findLastParticleInChain(static_cast<const reco::GenParticle*>(&*wd));
             addRadiation(motherPartIdx_, last_q, target);
             if (std::abs(wd->pdgId()) == TopDecayID::tauID) {
               // add tau daughters
               // currently it adds k-mesons etc as well, which
               // is not what we want.
               addDaughters(motherPartIdx_, wd, target);
             }
           }
         }
  
       } else {
         if (addRadiation_ && (td->pdgId() == TopDecayID::glueID || std::abs(td->pdgId()) < TopDecayID::bID)) {
           // collect additional radiation from the top
           std::unique_ptr<reco::GenParticle> radPtr(
               new reco::GenParticle(td->threeCharge(), td->p4(), td->vertex(), td->pdgId(), td->status(), false));
           target.push_back(*radPtr);
         }
         //other top daughters like Zq for FCNC
         // for pythia 6 many gluons end up here
         //std::cout << "other top daughters: to be implemented"
         //      << std::endl;
       }
     }
  
     // fill the map for the administration
     // of daughter indices
     refs_[iTop] = topDaughters;
     refs_[iW] = wDaughters;
   }
 }

References funct::abs(), addDaughters(), addRadiation(), addRadiation_, reco::Candidate::begin(), TopDecayID::bID, reco::Candidate::end(), fillMode_, findLastParticleInChain(), TopDecayID::glueID, kME, motherPartIdx_, refs_, OrderedSet::t, filterCSVwithJSON::target, TopDecayID::tauID, and TopDecayID::WID.

◆ fillListing() [2/2]

void TopDecaySubset::fillListing	(	const std::vector< const reco::GenParticle * > &	tops,
		reco::GenParticleCollection &	target
	)

private

fill output vector for full decay chain

Definition at line 308 of file TopDecaySubset.cc.

                                                                     {
   unsigned int statusFlag;
   // determine status flag of the new
   // particle depending on the FillMode
   fillMode_ == kME ? statusFlag = 3 : statusFlag = 2;
  
   for (std::vector<const reco::GenParticle*>::const_iterator it = tops.begin(); it != tops.end(); ++it) {
     const reco::GenParticle* t = *it;
     // if particle is top or anti-top
     std::unique_ptr<reco::GenParticle> topPtr(
         new reco::GenParticle(t->threeCharge(), p4(it, statusFlag), t->vertex(), t->pdgId(), statusFlag, false));
     target.push_back(*topPtr);
     ++motherPartIdx_;
     // keep the top index for the map to manage the daughter refs
     int iTop = motherPartIdx_;
     std::vector<int> topDaughters;
     // define the W boson index (to be set later) for the map to
     // manage the daughter refs
     int iW = 0;
     std::vector<int> wDaughters;
     // sanity check
     if (showerModel_ != kPythia && t->begin() == t->end())
       throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && t->begin()==t->end()\n");
     //iterate over top daughters
     for (reco::GenParticle::const_iterator td = ((showerModel_ == kPythia) ? t->begin() : t->begin()->begin());
          td != ((showerModel_ == kPythia) ? t->end() : t->begin()->end());
          ++td) {
       if (td->status() == TopDecayID::unfrag && std::abs(td->pdgId()) <= TopDecayID::bID) {
         // if particle is beauty or other quark
         std::unique_ptr<reco::GenParticle> bPtr(
             new reco::GenParticle(td->threeCharge(), p4(td, statusFlag), td->vertex(), td->pdgId(), statusFlag, false));
         target.push_back(*bPtr);
         // increment & push index of the top daughter
         topDaughters.push_back(++motherPartIdx_);
         if (addRadiation_) {
           addRadiation(motherPartIdx_, td, target);
         }
       }
       // sanity check
       if (showerModel_ != kPythia && td->begin() == td->end())
         throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && td->begin()==td->end()\n");
       reco::GenParticle::const_iterator buffer = (showerModel_ == kPythia) ? td : td->begin();
       if (buffer->status() == TopDecayID::unfrag && std::abs(buffer->pdgId()) == TopDecayID::WID) {
         // if particle is a W boson
         std::unique_ptr<reco::GenParticle> wPtr(new reco::GenParticle(
             buffer->threeCharge(), p4(buffer, statusFlag), buffer->vertex(), buffer->pdgId(), statusFlag, true));
         target.push_back(*wPtr);
         // increment & push index of the top daughter
         topDaughters.push_back(++motherPartIdx_);
         // keep the W idx for the map
         iW = motherPartIdx_;
         if (addRadiation_) {
           addRadiation(motherPartIdx_, buffer, target);
         }
         if (showerModel_ != kPythia && buffer->begin() == buffer->end())
           throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && buffer->begin()==buffer->end()\n");
         // iterate over W daughters
         for (reco::GenParticle::const_iterator wd =
                  ((showerModel_ == kPythia) ? buffer->begin() : buffer->begin()->begin());
              wd != ((showerModel_ == kPythia) ? buffer->end() : buffer->begin()->end());
              ++wd) {
           // make sure the W daughter is of status unfrag and not the W itself
           if (wd->status() == TopDecayID::unfrag && !(std::abs(wd->pdgId()) == TopDecayID::WID)) {
             std::unique_ptr<reco::GenParticle> qPtr(new reco::GenParticle(
                 wd->threeCharge(), p4(wd, statusFlag), wd->vertex(), wd->pdgId(), statusFlag, false));
             target.push_back(*qPtr);
             // increment & push index of the top daughter
             wDaughters.push_back(++motherPartIdx_);
             if (wd->status() == TopDecayID::unfrag && std::abs(wd->pdgId()) == TopDecayID::tauID) {
               // add tau daughters if the particle is a tau pass
               // the daughter of the tau which is of status 2
               //addDaughters(motherPartIdx_, wd->begin(), target);
               // add tau daughters if the particle is a tau pass
               // the tau itself, which may add a tau daughter of
               // of status 2 to the listing
               addDaughters(motherPartIdx_, wd, target);
             }
           }
         }
       }
       if (addRadiation_ && buffer->status() == TopDecayID::stable &&
           (buffer->pdgId() == TopDecayID::glueID || std::abs(buffer->pdgId()) < TopDecayID::bID)) {
         // collect additional radiation from the top
         std::unique_ptr<reco::GenParticle> radPtr(new reco::GenParticle(
             buffer->threeCharge(), buffer->p4(), buffer->vertex(), buffer->pdgId(), statusFlag, false));
         target.push_back(*radPtr);
         // increment & push index of the top daughter
         topDaughters.push_back(++motherPartIdx_);
       }
     }
     // add potential sisters of the top quark;
     // only for top to prevent double counting
     if (t->numberOfMothers() > 0 && t->pdgId() == TopDecayID::tID) {
       for (reco::GenParticle::const_iterator ts = t->mother()->begin(); ts != t->mother()->end(); ++ts) {
         // loop over all daughters of the top mother i.e.
         // the two top quarks and their potential sisters
         if (std::abs(ts->pdgId()) != t->pdgId() && ts->pdgId() != t->mother()->pdgId()) {
           // add all further particles but the two top quarks and potential
           // cases where the mother of the top has itself as daughter
           reco::GenParticle* cand =
               new reco::GenParticle(ts->threeCharge(), ts->p4(), ts->vertex(), ts->pdgId(), ts->status(), false);
           std::unique_ptr<reco::GenParticle> sPtr(cand);
           target.push_back(*sPtr);
           if (ts->begin() != ts->end()) {
             // in case the sister has daughters increment
             // and add the first generation of daughters
             addDaughters(++motherPartIdx_, ts->begin(), target, false);
           }
         }
       }
     }
     // fill the map for the administration
     // of daughter indices
     refs_[iTop] = topDaughters;
     refs_[iW] = wDaughters;
   }
 }

References funct::abs(), addDaughters(), addRadiation(), addRadiation_, TopDecayID::bID, edmScanValgrind::buffer, fillMode_, GenParticle::GenParticle, TopDecayID::glueID, kME, kPythia, edm::errors::LogicError, motherPartIdx_, p4(), refs_, showerModel_, TopDecayID::stable, OrderedSet::t, filterCSVwithJSON::target, TopDecayID::tauID, TopDecayID::tID, TopDecayID::unfrag, and TopDecayID::WID.

Referenced by produce().

◆ fillReferences()

void TopDecaySubset::fillReferences	(	const reco::GenParticleRefProd &	refProd,
		reco::GenParticleCollection &	target
	)

private

fill references for output vector

Definition at line 663 of file TopDecaySubset.cc.

                                                                                                      {
   int idx = 0;
   for (reco::GenParticleCollection::iterator p = sel.begin(); p != sel.end(); ++p, ++idx) {
     //find daughter reference vectors in refs_ and add daughters
     std::map<int, std::vector<int> >::const_iterator daughters = refs_.find(idx);
     if (daughters != refs_.end()) {
       for (std::vector<int>::const_iterator daughter = daughters->second.begin(); daughter != daughters->second.end();
            ++daughter) {
         reco::GenParticle* part = dynamic_cast<reco::GenParticle*>(&(*p));
         if (part == nullptr) {
           throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
         }
         part->addDaughter(reco::GenParticleRef(ref, *daughter));
         sel[*daughter].addMother(reco::GenParticleRef(ref, idx));
       }
     }
   }
 }

References BPhysicsValidation_cfi::daughters, Exception, training_settings::idx, edm::errors::InvalidReference, AlCaHLTBitMon_ParallelJobs::p, refs_, and EgammaValidation_Wenu_cff::sel.

Referenced by produce().

◆ findDecayingTops()

std::vector< const reco::GenParticle * > TopDecaySubset::findDecayingTops ( const reco::GenParticleCollection & parts )

private

find decaying top quarks (quarks that decay to qW) for Pythia6 this is identical to findPrimalTops

Definition at line 132 of file TopDecaySubset.cc.

                                                                                                          {
   std::vector<const reco::GenParticle*> tops;
   for (reco::GenParticleCollection::const_iterator t = parts.begin(); t != parts.end(); ++t) {
     if (std::abs(t->pdgId()) != TopDecayID::tID)
       continue;
  
     bool hasTopDaughter = false;
     for (unsigned idx = 0; idx < t->numberOfDaughters(); ++idx) {
       if (std::abs(t->daughter(idx)->pdgId()) == TopDecayID::tID)
         hasTopDaughter = true;
     }
  
     if (hasTopDaughter)  // not a decaying top
       continue;
     tops.push_back(&(*t));
   }
  
   return tops;
 }

References funct::abs(), training_settings::idx, contentValuesFiles::parts, OrderedSet::t, and TopDecayID::tID.

Referenced by produce().

◆ findLastParticleInChain()

const reco::GenParticle * TopDecaySubset::findLastParticleInChain ( const reco::GenParticle * p )

private

find W bosons that come from top quark decays and decay themselves (end of the MC chain) for Pythia6 this is identical to findPrimalW find the last particle in a (potentially) long chain of state transitions e.g. top[status==22]-> top[status==44 -> top[status==44] -> top[status==44] -> top[status==62] this function would pick the top with status 62

Definition at line 177 of file TopDecaySubset.cc.

                                                                                        {
   int particleID = std::abs(p->pdgId());
   bool containsItself = false;
   unsigned int d_idx = 0;
   for (unsigned idx = 0; idx < p->numberOfDaughters(); ++idx) {
     if (std::abs(p->daughter(idx)->pdgId()) == particleID) {
       containsItself = true;
       d_idx = idx;
     }
   }
  
   if (!containsItself)
     return p;
   else {
     if (showerModel_ == kPythia) {
       // Pythia6 has a weird idea of W bosons (and maybe other particles)
       // W (status == 3) -> q qbar' W. The new W is status 2 and has no daughters
       if (p->status() == 3)
         return p;
     }
     return findLastParticleInChain(static_cast<const reco::GenParticle*>(p->daughter(d_idx)));
   }
 }

References funct::abs(), training_settings::idx, kPythia, AlCaHLTBitMon_ParallelJobs::p, EgammaObjectsElectrons_cfi::particleID, and showerModel_.

Referenced by fillListing().

◆ findPrimalTops()

std::vector< const reco::GenParticle * > TopDecaySubset::findPrimalTops ( const reco::GenParticleCollection & parts )

private

find primal top quarks (top quarks from the hard interaction) for Pythia6 this is identical to findDecayingTops

Definition at line 110 of file TopDecaySubset.cc.

                                                                                                        {
   std::vector<const reco::GenParticle*> tops;
   for (reco::GenParticleCollection::const_iterator t = parts.begin(); t != parts.end(); ++t) {
     if (std::abs(t->pdgId()) != TopDecayID::tID)
       continue;
  
     bool hasTopMother = false;
     for (unsigned idx = 0; idx < t->numberOfMothers(); ++idx) {
       if (std::abs(t->mother(idx)->pdgId()) == TopDecayID::tID)
         hasTopMother = true;
     }
  
     if (hasTopMother)  // not a primal top
       continue;
     tops.push_back(&(*t));
   }
  
   return tops;
 }

References funct::abs(), training_settings::idx, contentValuesFiles::parts, OrderedSet::t, and TopDecayID::tID.

Referenced by produce().

◆ findPrimalW()

const reco::GenParticle * TopDecaySubset::findPrimalW ( const reco::GenParticle * top )

private

find W bosons that come from top quark decays for Pythia6 this is identical to findDecayingW

Definition at line 154 of file TopDecaySubset.cc.

                                                                              {
   unsigned int w_index = 0;
   for (unsigned idx = 0; idx < top->numberOfDaughters(); ++idx) {
     if (std::abs(top->daughter(idx)->pdgId()) == TopDecayID::WID) {
       w_index = idx;
       break;
     }
   }
   return static_cast<const reco::GenParticle*>(top->daughter(w_index));
 }

References funct::abs(), reco::CompositeRefCandidateT< D >::daughter(), training_settings::idx, reco::CompositeRefCandidateT< D >::numberOfDaughters(), reco::Candidate::pdgId(), and TopDecayID::WID.

◆ findTops()

std::vector< const reco::GenParticle * > TopDecaySubset::findTops ( const reco::GenParticleCollection & parts )

private

find top quarks in list of input particles

Definition at line 99 of file TopDecaySubset.cc.

                                                                                                  {
   std::vector<const reco::GenParticle*> tops;
   for (reco::GenParticleCollection::const_iterator t = parts.begin(); t != parts.end(); ++t) {
     if (std::abs(t->pdgId()) == TopDecayID::tID && t->status() == TopDecayID::unfrag)
       tops.push_back(&(*t));
   }
   return tops;
 }

References funct::abs(), contentValuesFiles::parts, OrderedSet::t, TopDecayID::tID, and TopDecayID::unfrag.

Referenced by produce().

◆ p4() [1/2]

reco::Particle::LorentzVector TopDecaySubset::p4	(	const reco::GenParticle::const_iterator	part,
		int	statusFlag
	)

private

calculate lorentz vector from input

calculate lorentz vector from input (dedicated to top reconstruction)

Definition at line 562 of file TopDecaySubset.cc.

                                                                                                      {
   // calculate the four vector for all top daughters from
   // their daughters including additional radiation
   if (statusFlag == TopDecayID::unfrag) {
     // return 4 momentum as it is
     return part->p4();
   }
   reco::Particle::LorentzVector vec;
   for (reco::GenParticle::const_iterator p = part->begin(); p != part->end(); ++p) {
     if (p->status() <= TopDecayID::stable && std::abs(p->pdgId()) == TopDecayID::WID) {
       vec = p->p4();
     } else {
       if (p->status() <= TopDecayID::stable) {
         // sum up the p4 of all stable particles
         // (of status 1 or 2)
         vec += p->p4();
       } else {
         if (p->status() == TopDecayID::unfrag) {
           // if the particle is unfragmented (i.e.
           // status 3) descend by one level
           vec += p4(p, statusFlag);
         }
       }
     }
   }
   return vec;
 }

References funct::abs(), AlCaHLTBitMon_ParallelJobs::p, p4(), TopDecayID::stable, TopDecayID::unfrag, and TopDecayID::WID.

Referenced by Tau.Tau::dxy_approx(), Tau.Tau::dz(), and Lepton.Lepton::p4WithFSR().

◆ p4() [2/2]

reco::Particle::LorentzVector TopDecaySubset::p4	(	const std::vector< const reco::GenParticle * >::const_iterator	top,
		int	statusFlag
	)

private

calculate lorentz vector from input (dedicated to top reconstruction)

calculate lorentz vector from input

Definition at line 524 of file TopDecaySubset.cc.

                                                                  {
   // calculate the four vector for top/anti-top quarks from
   // the W boson and the b quark plain or including all
   // additional radiation depending on switch 'plain'
   if (statusFlag == TopDecayID::unfrag) {
     // return 4 momentum as it is
     return (*top)->p4();
   }
   reco::Particle::LorentzVector vec;
   for (reco::GenParticle::const_iterator p = (*top)->begin(); p != (*top)->end(); ++p) {
     if (p->status() == TopDecayID::unfrag) {
       // descend by one level for each
       // status 3 particle on the way
       vec += p4(p, statusFlag);
     } else {
       if (std::abs((*top)->pdgId()) == TopDecayID::WID) {
         // in case of a W boson skip the status
         // 2 particle to prevent double counting
         if (std::abs(p->pdgId()) != TopDecayID::WID)
           vec += p->p4();
       } else {
         // add all four vectors for each stable
         // particle (status 1 or 2) on the way
         vec += p->p4();
         if (vec.mass() - (*top)->mass() > 0) {
           // continue adding up gluons and qqbar pairs on the top
           // line untill the nominal top mass is reached; then
           // break in order to prevent picking up virtualities
           break;
         }
       }
     }
   }
   return vec;
 }

References funct::abs(), AlCaHLTBitMon_ParallelJobs::p, TopDecayID::unfrag, and TopDecayID::WID.

Referenced by Tau.Tau::dxy_approx(), Tau.Tau::dz(), fillListing(), p4(), and Lepton.Lepton::p4WithFSR().

◆ produce()

void TopDecaySubset::produce	(	edm::Event &	event,
		const edm::EventSetup &	setup
	)

overridevirtual

write output into the event

Implements edm::EDProducer.

Definition at line 44 of file TopDecaySubset.cc.

                                                                         {
   // create target vector
   std::unique_ptr<reco::GenParticleCollection> target(new reco::GenParticleCollection);
  
   // get source collection
   edm::Handle<reco::GenParticleCollection> src;
   event.getByToken(srcToken_, src);
  
   // find out what generator we are dealing with
   if (showerModel_ == kStart && runMode_ == kRun2) {
     showerModel_ = checkShowerModel(event);
   }
  
   // find top quarks in list of input particles
   std::vector<const reco::GenParticle*> tops;
   if (runMode_ == kRun1)
     tops = findTops(*src);
   else
     tops = findPrimalTops(*src);
  
   // determine shower model (only in first event)
   if (showerModel_ == kStart && runMode_ == kRun1)
     showerModel_ = checkShowerModel(tops);
  
   if (showerModel_ != kNone) {
     // check sanity of W bosons
     if (runMode_ == kRun1)
       checkWBosons(tops);
     else {
       // nothing for the moment
     }
  
     // get ref product from the event
     const reco::GenParticleRefProd ref = event.getRefBeforePut<reco::GenParticleCollection>();
     // clear existing refs
     clearReferences();
     if (runMode_ == kRun1) {
       // fill output
       fillListing(tops, *target);
       // fill references
       fillReferences(ref, *target);
     } else {
       std::vector<const reco::GenParticle*> decaying_tops = findDecayingTops(*src);
       // fill output
       fillListing(tops, decaying_tops, *target);
       // fill references
       fillReferences(ref, *target);
     }
   }
  
   // write vectors to the event
   event.put(std::move(target));
 }

References checkShowerModel(), checkWBosons(), clearReferences(), fillListing(), fillReferences(), findDecayingTops(), findPrimalTops(), findTops(), kNone, kRun1, kRun2, kStart, eostools::move(), runMode_, showerModel_, TrackRefitter_38T_cff::src, srcToken_, and filterCSVwithJSON::target.

Member Data Documentation

◆ addRadiation_

bool TopDecaySubset::addRadiation_

private

add radiation or not?

Definition at line 100 of file TopDecaySubset.h.

Referenced by fillListing().

◆ fillMode_

FillMode TopDecaySubset::fillMode_

private

print the whole list of input particles or not? mode of decaySubset creation

Definition at line 103 of file TopDecaySubset.h.

Referenced by fillListing(), and TopDecaySubset().

◆ genEventInfo_srcToken_

edm::EDGetTokenT<GenEventInfoProduct> TopDecaySubset::genEventInfo_srcToken_

private

input tag for the genEventInfo source

Definition at line 98 of file TopDecaySubset.h.

Referenced by checkShowerModel().

◆ motherPartIdx_

int TopDecaySubset::motherPartIdx_

private

index in new evt listing of parts with daughters; has to be set to -1 in produce to deliver consistent results!

Definition at line 112 of file TopDecaySubset.h.

Referenced by clearReferences(), and fillListing().

◆ refs_

std::map<int, std::vector<int> > TopDecaySubset::refs_

private

management of daughter indices for fillRefs

Definition at line 114 of file TopDecaySubset.h.

Referenced by addDaughters(), addRadiation(), clearReferences(), fillListing(), and fillReferences().

◆ runMode_

RunMode TopDecaySubset::runMode_

private

run mode (Run1 || Run2)

Definition at line 107 of file TopDecaySubset.h.

Referenced by produce(), and TopDecaySubset().

◆ showerModel_

ShowerModel TopDecaySubset::showerModel_

private

parton shower mode (filled in checkShowerModel)

Definition at line 105 of file TopDecaySubset.h.

Referenced by checkWBosons(), fillListing(), findLastParticleInChain(), and produce().

◆ srcToken_

edm::EDGetTokenT<reco::GenParticleCollection> TopDecaySubset::srcToken_

private

input tag for the genParticle source

Definition at line 96 of file TopDecaySubset.h.

Referenced by produce().

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Enumeration Documentation

◆ FillMode

◆ RunMode

◆ ShowerModel

Constructor & Destructor Documentation

◆ TopDecaySubset()

◆ ~TopDecaySubset()

Member Function Documentation

◆ addDaughters()

◆ addRadiation() [1/2]

◆ addRadiation() [2/2]

◆ checkShowerModel() [1/2]

◆ checkShowerModel() [2/2]

◆ checkWBosons()

◆ clearReferences()

◆ fillListing() [1/2]

◆ fillListing() [2/2]

◆ fillReferences()

◆ findDecayingTops()

◆ findLastParticleInChain()

◆ findPrimalTops()

◆ findPrimalW()

◆ findTops()

◆ p4() [1/2]

◆ p4() [2/2]

◆ produce()

Member Data Documentation

◆ addRadiation_

◆ fillMode_

◆ genEventInfo_srcToken_

◆ motherPartIdx_

◆ refs_

◆ runMode_

◆ showerModel_

◆ srcToken_