Clusters hadrons, partons, and jet contituents to determine the jet flavour. More...

#include <PhysicsTools/JetMCAlgos/plugins/JetFlavourClustering.cc>

Inheritance diagram for JetFlavourClustering:

Public Member Functions
	JetFlavourClustering (const edm::ParameterSet &)

	~JetFlavourClustering ()

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)

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

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

Static Public Member Functions inherited from edm::stream::EDProducerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Private Member Functions
void	assignToSubjets (const reco::GenParticleRefVector &clusteredParticles, const edm::Handle< edm::View< reco::Jet > > &subjets, const std::vector< int > &subjetIndices, std::vector< reco::GenParticleRefVector > &assignedParticles)

void	insertGhosts (const edm::Handle< reco::GenParticleRefVector > &particles, const double ghostRescaling, const bool isHadron, const bool isbHadron, const bool isParton, const bool isLepton, std::vector< fastjet::PseudoJet > &constituents)

void	matchGroomedJets (const edm::Handle< edm::View< reco::Jet > > &jets, const edm::Handle< edm::View< reco::Jet > > &matchedJets, std::vector< int > &matchedIndices)

void	matchReclusteredJets (const edm::Handle< edm::View< reco::Jet > > &jets, const std::vector< fastjet::PseudoJet > &matchedJets, std::vector< int > &matchedIndices)

void	matchSubjets (const std::vector< int > &groomedIndices, const edm::Handle< edm::View< reco::Jet > > &groomedJets, const edm::Handle< edm::View< reco::Jet > > &subjets, std::vector< std::vector< int > > &matchedIndices)

virtual void	produce (edm::Event &, const edm::EventSetup &)

void	setFlavours (const reco::GenParticleRefVector &clusteredbHadrons, const reco::GenParticleRefVector &clusteredcHadrons, const reco::GenParticleRefVector &clusteredPartons, int &hadronFlavour, int &partonFlavour)

Private Attributes
const edm::EDGetTokenT < reco::GenParticleRefVector >	bHadronsToken_

const edm::EDGetTokenT < reco::GenParticleRefVector >	cHadronsToken_

ClusterSequencePtr	fjClusterSeq_

JetDefPtr	fjJetDefinition_

const double	ghostRescaling_

edm::EDGetTokenT< edm::View < reco::Jet > >	groomedJetsToken_

const bool	hadronFlavourHasPriority_

const std::string	jetAlgorithm_

const double	jetPtMin_

const edm::EDGetTokenT < edm::View< reco::Jet > >	jetsToken_

edm::EDGetTokenT < reco::GenParticleRefVector >	leptonsToken_

const edm::EDGetTokenT < reco::GenParticleRefVector >	partonsToken_

const double	relPtTolerance_

const double	rParam_

edm::EDGetTokenT< edm::View < reco::Jet > >	subjetsToken_

const bool	useLeptons_

const bool	useSubjets_

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

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

Clusters hadrons, partons, and jet contituents to determine the jet flavour.

This producer clusters hadrons, partons and jet contituents to determine the jet flavour. The jet flavour information is stored in the event as an AssociationVector which associates an object of type JetFlavourInfo to each of the jets.

The producer takes as input jets and hadron and partons selected by the HadronAndPartonSelector producer. The hadron and parton four-momenta are rescaled by a very small number (default rescale factor is 10e-18) which turns them into the so-called "ghosts". The "ghost" hadrons and partons are clustered together with all of the jet constituents. It is important to use the same clustering algorithm and jet size as for the original input jet collection. Since the "ghost" hadrons and partons are extremely soft, the resulting jet collection will be practically identical to the original one but now with "ghost" hadrons and partons clustered inside jets. The jet flavour is determined based on the "ghost" hadrons clustered inside a jet:

jet is considered a b jet if there is at least one b "ghost" hadron clustered inside it (hadronFlavour = 5)
jet is considered a c jet if there is at least one c and no b "ghost" hadrons clustered inside it (hadronFlavour = 4)
jet is considered a light-flavour jet if there are no b or c "ghost" hadrons clustered inside it (hadronFlavour = 0)

To further assign a more specific flavour to light-flavour jets, "ghost" partons are used:

jet is considered a b jet if there is at least one b "ghost" parton clustered inside it (partonFlavour = 5)
jet is considered a c jet if there is at least one c and no b "ghost" partons clustered inside it (partonFlavour = 4)
jet is considered a light-flavour jet if there are light-flavour and no b or c "ghost" partons clustered inside it. The jet is assigned the flavour of the hardest light-flavour "ghost" parton clustered inside it (partonFlavour = 1, 2, 3, or 21)
jet has an undefined flavour if there are no "ghost" partons clustered inside it (partonFlavour = 0)

In rare instances a conflict between the hadron- and parton-based flavours can occur. In such cases it is possible to keep both flavours or to give priority to the hadron-based flavour. This is controlled by the 'hadronFlavourHasPriority' switch. The priority is given to the hadron-based flavour as follows:

if hadronFlavour==0 && (partonFlavour==4 || partonFlavour==5): partonFlavour is set to the flavour of the hardest light-flavour parton clustered inside the jet if such parton exists. Otherwise, the parton flavour is left undefined
if hadronFlavour!=0 && hadronFlavour!=partonFlavour: partonFlavour is set equal to hadronFlavour

The producer is also capable of assigning the flavour to subjets of fat jets, in which case it produces an additional AssociationVector providing the flavour information for subjets. In order to assign the flavour to subjets, three input jet collections are required:

jets, in this case represented by fat jets
groomed jets, which is a collection of fat jets from which the subjets are derived (e.g. pruned, filtered, soft drop, top-tagged, etc. jets)
subjets, derived from the groomed fat jets

The "ghost" hadrons and partons clustered inside a fat jet are assigned to the closest subjet in the rapidity-phi space. Once hadrons and partons have been assigned to subjets, the subjet flavour is determined in the same way as for jets. The reason for requiring three jet collections as input in order to determine the subjet flavour is to avoid possible inconsistencies between the fat jet and subjet flavours (such as a non-b fat jet having a b subjet and vice versa) as well as the fact that re-clustering the constituents of groomed fat jets will generally result in a jet collection different from the input groomed fat jets. Also note that "ghost" particles generally cannot be clustered inside subjets in the same way this is done for fat jets. This is because some of the jet grooming techniques could reject such very soft particle. So instead, the "ghost" particles are assigned to the closest subjet.

Finally, "ghost" leptons can also be clustered inside jets but they are not used in any way to determine the jet flavour. This functionality is optional and is potentially useful to identify jets from hadronic taus.

For more details, please refer to https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideBTagMCTools

Definition at line 141 of file JetFlavourClustering.cc.

Constructor & Destructor Documentation

JetFlavourClustering::JetFlavourClustering ( const edm::ParameterSet & iConfig )

explicit

Definition at line 207 of file JetFlavourClustering.cc.

References Exception, fjJetDefinition_, edm::ParameterSet::getParameter(), groomedJetsToken_, jetAlgorithm_, leptonsToken_, rParam_, subjetsToken_, useLeptons_, and useSubjets_.

                                                                          :
 
    jetsToken_(consumes<edm::View<reco::Jet> >( iConfig.getParameter<edm::InputTag>("jets")) ),
    bHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("bHadrons") )),
    cHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("cHadrons") )),
    partonsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("partons") )),
    jetAlgorithm_(iConfig.getParameter<std::string>("jetAlgorithm")),
    rParam_(iConfig.getParameter<double>("rParam")),
    jetPtMin_(0.), // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
    ghostRescaling_(iConfig.exists("ghostRescaling") ? iConfig.getParameter<double>("ghostRescaling") : 1e-18),
    relPtTolerance_(iConfig.exists("relPtTolerance") ? iConfig.getParameter<double>("relPtTolerance") : 1e-03), // 0.1% relative difference in Pt should be sufficient to detect possible misconfigurations
    hadronFlavourHasPriority_(iConfig.getParameter<bool>("hadronFlavourHasPriority")),
    useSubjets_(iConfig.exists("groomedJets") && iConfig.exists("subjets")),
    useLeptons_(iConfig.exists("leptons"))
 
 {
    // register your products
    produces<reco::JetFlavourInfoMatchingCollection>();
    if( useSubjets_ )
      produces<reco::JetFlavourInfoMatchingCollection>("SubJets");
 
    // set jet algorithm
    if (jetAlgorithm_=="Kt")
      fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::kt_algorithm, rParam_) );
    else if (jetAlgorithm_=="CambridgeAachen")
      fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam_) );
    else if (jetAlgorithm_=="AntiKt")
      fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam_) );
    else
      throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm_ << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
 
    if (useSubjets_) {
      groomedJetsToken_=consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("groomedJets"));
      subjetsToken_=consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("subjets"));
    }
    if ( useLeptons_ ) {
      leptonsToken_=consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("leptons") );
    }
 
 }

JetFlavourClustering::~JetFlavourClustering ( )

Definition at line 249 of file JetFlavourClustering.cc.

 {
 
    // do anything here that needs to be done at desctruction time
    // (e.g. close files, deallocate resources etc.)
 
 }

Member Function Documentation

void JetFlavourClustering::assignToSubjets	(	const reco::GenParticleRefVector &	clusteredParticles,
		const edm::Handle< edm::View< reco::Jet > > &	subjets,
		const std::vector< int > &	subjetIndices,
		std::vector< reco::GenParticleRefVector > &	assignedParticles
	)

private

Definition at line 697 of file JetFlavourClustering.cc.

References edm::RefVector< C, T, F >::begin(), reco::deltaR2(), HLT_FULL_cff::distance, edm::RefVector< C, T, F >::end(), and phi.

Referenced by produce().

 {
    // loop over clustered particles and assign them to different subjets based on smallest dR
    for(reco::GenParticleRefVector::const_iterator it = clusteredParticles.begin(); it != clusteredParticles.end(); ++it)
    {
      std::vector<double> dR2toSubjets;
 
      for(size_t sj=0; sj<subjetIndices.size(); ++sj)
        dR2toSubjets.push_back( reco::deltaR2( (*it)->rapidity(), (*it)->phi(), subjets->at(subjetIndices.at(sj)).rapidity(), subjets->at(subjetIndices.at(sj)).phi() ) );
 
      // find the closest subjet
      int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
 
      assignedParticles.at(closestSubjetIdx).push_back( *it );
    }
 }

void JetFlavourClustering::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

static

Definition at line 719 of file JetFlavourClustering.cc.

References edm::ConfigurationDescriptions::addDefault(), and edm::ParameterSetDescription::setUnknown().

                                                                                  {
   //The following says we do not know what parameters are allowed so do no validation
   // Please change this to state exactly what you do use, even if it is no parameters
   edm::ParameterSetDescription desc;
   desc.setUnknown();
   descriptions.addDefault(desc);
 }

void JetFlavourClustering::insertGhosts	(	const edm::Handle< reco::GenParticleRefVector > &	particles,
		const double	ghostRescaling,
		const bool	isHadron,
		const bool	isbHadron,
		const bool	isParton,
		const bool	isLepton,
		std::vector< fastjet::PseudoJet > &	constituents
	)

private

Definition at line 485 of file JetFlavourClustering.cc.

References AlCaHLTBitMon_ParallelJobs::p.

Referenced by produce().

 {
    // insert "ghost" particles in the vector of jet constituents
    for(reco::GenParticleRefVector::const_iterator it = particles->begin(); it != particles->end(); ++it)
    {
      if((*it)->pt() == 0)
      {
        edm::LogInfo("NullTransverseMomentum") << "dropping input ghost candidate with pt=0";
        continue;
      }
      fastjet::PseudoJet p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
      p*=ghostRescaling; // rescale particle momentum
      p.set_user_info(new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
      constituents.push_back(p);
    }
 }

void JetFlavourClustering::matchGroomedJets	(	const edm::Handle< edm::View< reco::Jet > > &	jets,
		const edm::Handle< edm::View< reco::Jet > > &	matchedJets,
		std::vector< int > &	matchedIndices
	)

private

Definition at line 549 of file JetFlavourClustering.cc.

References reco::deltaR2(), HLT_FULL_cff::distance, spr::find(), j, fwrapper::jets, rParam_, and mathSSE::sqrt().

Referenced by produce().

 {
    std::vector<bool> jetLocks(jets->size(),false);
    std::vector<int>  jetIndices;
 
    for(size_t gj=0; gj<groomedJets->size(); ++gj)
    {
      double matchedDR2 = 1e9;
      int matchedIdx = -1;
 
      if( groomedJets->at(gj).pt()>0. ) // skip pathological cases of groomed jets with Pt=0
      {
        for(size_t j=0; j<jets->size(); ++j)
        {
          if( jetLocks.at(j) ) continue; // skip jets that have already been matched
 
          double tempDR2 = reco::deltaR2( jets->at(j).rapidity(), jets->at(j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi() );
          if( tempDR2 < matchedDR2 )
          {
            matchedDR2 = tempDR2;
            matchedIdx = j;
          }
        }
      }
 
      if( matchedIdx>=0 )
      {
        if ( matchedDR2 > rParam_*rParam_ )
        {
          edm::LogWarning("MatchedJetsFarApart") << "Matched groomed jet " << gj << " and original jet " << matchedIdx <<" are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam_ << ".\n"
                                                 << "This is not expected so the matching of these two jets has been discarded. Please check that the two jet collections belong to each other.";
          matchedIdx = -1;
        }
        else
          jetLocks.at(matchedIdx) = true;
      }
      jetIndices.push_back(matchedIdx);
    }
 
    for(size_t j=0; j<jets->size(); ++j)
    {
      std::vector<int>::iterator matchedIndex = std::find( jetIndices.begin(), jetIndices.end(), j );
 
      matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
    }
 }

void JetFlavourClustering::matchReclusteredJets	(	const edm::Handle< edm::View< reco::Jet > > &	jets,
		const std::vector< fastjet::PseudoJet > &	matchedJets,
		std::vector< int > &	matchedIndices
	)

private

Definition at line 507 of file JetFlavourClustering.cc.

References reco::deltaR2(), j, fwrapper::jets, rParam_, and mathSSE::sqrt().

Referenced by produce().

 {
    std::vector<bool> matchedLocks(reclusteredJets.size(),false);
 
    for(size_t j=0; j<jets->size(); ++j)
    {
      double matchedDR2 = 1e9;
      int matchedIdx = -1;
 
      for(size_t rj=0; rj<reclusteredJets.size(); ++rj)
      {
        if( matchedLocks.at(rj) ) continue; // skip jets that have already been matched
 
        double tempDR2 = reco::deltaR2( jets->at(j).rapidity(), jets->at(j).phi(), reclusteredJets.at(rj).rapidity(), reclusteredJets.at(rj).phi_std() );
        if( tempDR2 < matchedDR2 )
        {
          matchedDR2 = tempDR2;
          matchedIdx = rj;
        }
      }
 
      if( matchedIdx>=0 )
      {
        if ( matchedDR2 > rParam_*rParam_ )
        {
          edm::LogError("JetMatchingFailed") << "Matched reclustered jet " << matchedIdx << " and original jet " << j <<" are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam_ << ".\n"
                                             << "This is not expected so please check that the jet algorithm and jet size match those used for the original jet collection.";
        }
        else
          matchedLocks.at(matchedIdx) = true;
      }
      else
        edm::LogError("JetMatchingFailed") << "Matching reclustered to original jets failed. Please check that the jet algorithm and jet size match those used for the original jet collection.";
 
      matchedIndices.push_back(matchedIdx);
    }
 }

void JetFlavourClustering::matchSubjets	(	const std::vector< int > &	groomedIndices,
		const edm::Handle< edm::View< reco::Jet > > &	groomedJets,
		const edm::Handle< edm::View< reco::Jet > > &	subjets,
		std::vector< std::vector< int > > &	matchedIndices
	)

private

Definition at line 600 of file JetFlavourClustering.cc.

References g, and alignCSCRings::s.

Referenced by produce().

 {
    for(size_t g=0; g<groomedIndices.size(); ++g)
    {
      std::vector<int> subjetIndices;
 
      if( groomedIndices.at(g)>=0 )
      {
        for(size_t s=0; s<groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s)
        {
          const edm::Ptr<reco::Candidate> & subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
 
          for(size_t sj=0; sj<subjets->size(); ++sj)
          {
            if( subjet == edm::Ptr<reco::Candidate>(subjets->ptrAt(sj)) )
            {
              subjetIndices.push_back(sj);
              break;
            }
          }
        }
 
        if( subjetIndices.size() == 0 )
          edm::LogError("SubjetMatchingFailed") << "Matching subjets to original jets failed. Please check that the groomed jet and subjet collections belong to each other.";
 
        matchedIndices.push_back(subjetIndices);
      }
      else
        matchedIndices.push_back(subjetIndices);
    }
 }

void JetFlavourClustering::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

privatevirtual

Implements edm::stream::EDProducerBase.

Definition at line 264 of file JetFlavourClustering.cc.

References funct::abs(), assignToSubjets(), bHadronsToken_, cHadronsToken_, fjClusterSeq_, fjJetDefinition_, edm::Event::getByToken(), ghostRescaling_, groomedJetsToken_, i, insertGhosts(), edm::Ptr< T >::isAvailable(), GhostInfo::isbHadron(), GhostInfo::isHadron(), GhostInfo::isLepton(), edm::Ptr< T >::isNonnull(), GhostInfo::isParton(), jetPtMin_, fwrapper::jets, jetsToken_, HLT_FULL_cff::leptons, leptonsToken_, visualization-live-secondInstance_cfg::m, matchGroomedJets(), matchReclusteredJets(), matchSubjets(), GhostInfo::particleRef(), partonsToken_, EnergyCorrector::pt, edm::RefVector< C, T, F >::push_back(), edm::Event::put(), relPtTolerance_, setFlavours(), subjetsToken_, useLeptons_, and useSubjets_.

 {
    edm::Handle<edm::View<reco::Jet> > jets;
    iEvent.getByToken(jetsToken_, jets);
 
    edm::Handle<edm::View<reco::Jet> > groomedJets;
    edm::Handle<edm::View<reco::Jet> > subjets;
    if( useSubjets_ )
    {
      iEvent.getByToken(groomedJetsToken_, groomedJets);
      iEvent.getByToken(subjetsToken_, subjets);
    }
 
    edm::Handle<reco::GenParticleRefVector> bHadrons;
    iEvent.getByToken(bHadronsToken_, bHadrons);
 
    edm::Handle<reco::GenParticleRefVector> cHadrons;
    iEvent.getByToken(cHadronsToken_, cHadrons);
 
    edm::Handle<reco::GenParticleRefVector> partons;
    iEvent.getByToken(partonsToken_, partons);
 
    edm::Handle<reco::GenParticleRefVector> leptons;
    if( useLeptons_ )
      iEvent.getByToken(leptonsToken_, leptons);
 
    std::auto_ptr<reco::JetFlavourInfoMatchingCollection> jetFlavourInfos( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(jets)) );
    std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
    if( useSubjets_ )
      subjetFlavourInfos = std::auto_ptr<reco::JetFlavourInfoMatchingCollection>( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(subjets)) );
 
    // vector of constituents for reclustering jets and "ghosts"
    std::vector<fastjet::PseudoJet> fjInputs;
    // loop over all input jets and collect all their constituents
    for(edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); ++it)
    {
      std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
      std::vector<edm::Ptr<reco::Candidate> >::const_iterator m;
      for( m = constituents.begin(); m != constituents.end(); ++m )
      {
        reco::CandidatePtr constit = *m;
        if(!constit.isNonnull() || !constit.isAvailable()) {
          edm::LogError("MissingJetConstituent") << "Jet constituent required for jet reclustering is missing. Reclustered jets are not guaranteed to reproduce the original jets!";
          continue;
        }
        if(constit->pt() == 0)
        {
          edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
          continue;
        }
        fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
      }
    }
    // insert "ghost" b hadrons in the vector of constituents
    insertGhosts(bHadrons, ghostRescaling_, true, true, false, false, fjInputs);
    // insert "ghost" c hadrons in the vector of constituents
    insertGhosts(cHadrons, ghostRescaling_, true, false, false, false, fjInputs);
    // insert "ghost" partons in the vector of constituents
    insertGhosts(partons, ghostRescaling_, false, false, true, false, fjInputs);
    // if used, insert "ghost" leptons in the vector of constituents
    if( useLeptons_ )
      insertGhosts(leptons, ghostRescaling_, false, false, false, true, fjInputs);
 
    // define jet clustering sequence
    fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs, *fjJetDefinition_ ) );
    // recluster jet constituents and inserted "ghosts"
    std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt( fjClusterSeq_->inclusive_jets(jetPtMin_) );
 
    if( inclusiveJets.size() < jets->size() )
      edm::LogError("TooFewReclusteredJets") << "There are fewer reclustered (" << inclusiveJets.size() << ") than original jets (" << jets->size() << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
 
    // match reclustered and original jets
    std::vector<int> reclusteredIndices;
    matchReclusteredJets(jets,inclusiveJets,reclusteredIndices);
 
    // match groomed and original jets
    std::vector<int> groomedIndices;
    if( useSubjets_ )
    {
      if( groomedJets->size() > jets->size() )
        edm::LogError("TooManyGroomedJets") << "There are more groomed (" << groomedJets->size() << ") than original jets (" << jets->size() << "). Please check that the two jet collections belong to each other.";
 
      matchGroomedJets(jets,groomedJets,groomedIndices);
    }
 
    // match subjets and original jets
    std::vector<std::vector<int> > subjetIndices;
    if( useSubjets_ )
    {
      matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
    }
 
    // determine jet flavour
    for(size_t i=0; i<jets->size(); ++i)
    {
      reco::GenParticleRefVector clusteredbHadrons;
      reco::GenParticleRefVector clusteredcHadrons;
      reco::GenParticleRefVector clusteredPartons;
      reco::GenParticleRefVector clusteredLeptons;
 
      // if matching reclustered to original jets failed
      if( reclusteredIndices.at(i) < 0 )
      {
        // set an empty JetFlavourInfo for this jet
        (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
      }
      else if( jets->at(i).pt() == 0 )
      {
        edm::LogWarning("NullTransverseMomentum") << "The original jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
 
        // set an empty JetFlavourInfo for this jet
        (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
 
        // if subjets are used
        if( useSubjets_ && subjetIndices.at(i).size()>0 )
        {
          // loop over subjets
          for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
          {
            // set an empty JetFlavourInfo for this subjet
            (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] = reco::JetFlavourInfo(reco::GenParticleRefVector(), reco::GenParticleRefVector(), reco::GenParticleRefVector(), reco::GenParticleRefVector(), 0, 0);
          }
        }
      }
      else
      {
        // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
        if( ( std::abs( inclusiveJets.at(reclusteredIndices.at(i)).pt() - jets->at(i).pt() ) / jets->at(i).pt() ) > relPtTolerance_ )
        {
          if( jets->at(i).pt() < 10. )  // special handling for low-Pt jets (Pt<10 GeV)
            edm::LogWarning("JetPtMismatchAtLowPt") << "The reclustered and original jet " << i << " have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << jets->at(i).pt() << " GeV, respectively).\n"
                                                    << "Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
                                                    << "Since the mismatch is at low Pt (Pt<10 GeV), it is ignored and only a warning is issued.\n"
                                                    << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
          else
            edm::LogError("JetPtMismatch") << "The reclustered and original jet " << i << " have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << jets->at(i).pt() << " GeV, respectively).\n"
                                           << "Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n"
                                           << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
        }
 
        // get jet constituents (sorted by Pt)
        std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(i)).constituents() );
 
        // loop over jet constituents and try to find "ghosts"
        for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
        {
          if( !it->has_user_info() ) continue; // skip if not a "ghost"
 
          // "ghost" hadron
          if( it->user_info<GhostInfo>().isHadron() )
          {
            // "ghost" b hadron
            if( it->user_info<GhostInfo>().isbHadron() )
              clusteredbHadrons.push_back(it->user_info<GhostInfo>().particleRef());
            // "ghost" c hadron
            else
              clusteredcHadrons.push_back(it->user_info<GhostInfo>().particleRef());
          }
          // "ghost" parton
          else if( it->user_info<GhostInfo>().isParton() )
            clusteredPartons.push_back(it->user_info<GhostInfo>().particleRef());
          // "ghost" lepton
          else if( it->user_info<GhostInfo>().isLepton() )
            clusteredLeptons.push_back(it->user_info<GhostInfo>().particleRef());
        }
 
        int hadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
        int partonFlavour = 0; // default parton flavour set to 0 (= undefined)
 
        // set hadron- and parton-based flavours
        setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
 
        // set the JetFlavourInfo for this jet
        (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
      }
 
      // if subjets are used, determine their flavour
      if( useSubjets_ )
      {
        if( subjetIndices.at(i).size()==0 ) continue; // continue if the original jet does not have subjets assigned
 
        // define vectors of GenParticleRefVectors for hadrons and partons assigned to different subjets
        std::vector<reco::GenParticleRefVector> assignedbHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
        std::vector<reco::GenParticleRefVector> assignedcHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
        std::vector<reco::GenParticleRefVector> assignedPartons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
        std::vector<reco::GenParticleRefVector> assignedLeptons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
 
        // loop over clustered b hadrons and assign them to different subjets based on smallest dR
        assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(i), assignedbHadrons);
        // loop over clustered c hadrons and assign them to different subjets based on smallest dR
        assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(i), assignedcHadrons);
        // loop over clustered partons and assign them to different subjets based on smallest dR
        assignToSubjets(clusteredPartons, subjets, subjetIndices.at(i), assignedPartons);
        // if used, loop over clustered leptons and assign them to different subjets based on smallest dR
        if( useLeptons_ )
          assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(i), assignedLeptons);
 
        // loop over subjets and determine their flavour
        for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
        {
          int subjetHadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
          int subjetPartonFlavour = 0; // default parton flavour set to 0 (= undefined)
 
          // set hadron- and parton-based flavours
          setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
 
          // set the JetFlavourInfo for this subjet
          (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] = reco::JetFlavourInfo(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), assignedLeptons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
        }
      }
    }
 
    // put jet flavour infos in the event
    iEvent.put( jetFlavourInfos );
    // put subjet flavour infos in the event
    if( useSubjets_ )
      iEvent.put( subjetFlavourInfos, "SubJets" );
 }

void JetFlavourClustering::setFlavours	(	const reco::GenParticleRefVector &	clusteredbHadrons,
		const reco::GenParticleRefVector &	clusteredcHadrons,
		const reco::GenParticleRefVector &	clusteredPartons,
		int &	hadronFlavour,
		int &	partonFlavour
	)

private

Definition at line 637 of file JetFlavourClustering.cc.

References funct::abs(), edm::RefVector< C, T, F >::begin(), edm::RefVector< C, T, F >::end(), hadronFlavourHasPriority_, CandMCTagUtils::isLightParton(), edm::Ref< C, T, F >::isNonnull(), edm::Ref< C, T, F >::isNull(), and edm::RefVector< C, T, F >::size().

Referenced by produce().

 {
    reco::GenParticleRef hardestParton;
    reco::GenParticleRef hardestLightParton;
    reco::GenParticleRef flavourParton;
 
    // loop over clustered partons (already sorted by Pt)
    for(reco::GenParticleRefVector::const_iterator it = clusteredPartons.begin(); it != clusteredPartons.end(); ++it)
    {
      // hardest parton
      if( hardestParton.isNull() )
        hardestParton = (*it);
      // hardest light-flavour parton
      if( hardestLightParton.isNull() )
      {
        if( CandMCTagUtils::isLightParton( *(*it) ) )
          hardestLightParton = (*it);
      }
      // c flavour
      if( flavourParton.isNull() && ( std::abs( (*it)->pdgId() ) == 4 ) )
        flavourParton = (*it);
      // b flavour gets priority
      if( std::abs( (*it)->pdgId() ) == 5 )
      {
        if( flavourParton.isNull() )
          flavourParton = (*it);
        else if( std::abs( flavourParton->pdgId() ) != 5 )
          flavourParton = (*it);
      }
    }
 
    // set hadron-based flavour
    if( clusteredbHadrons.size()>0 )
      hadronFlavour = 5;
    else if( clusteredcHadrons.size()>0 && clusteredbHadrons.size()==0 )
      hadronFlavour = 4;
    // set parton-based flavour
    if( flavourParton.isNull() )
    {
      if( hardestParton.isNonnull() ) partonFlavour = hardestParton->pdgId();
    }
    else
      partonFlavour = flavourParton->pdgId();
 
    // if enabled, check for conflicts between hadron- and parton-based flavours and give priority to the hadron-based flavour
    if( hadronFlavourHasPriority_ )
    {
      if( hadronFlavour==0 && (std::abs(partonFlavour)==4 || std::abs(partonFlavour)==5) )
        partonFlavour = ( hardestLightParton.isNonnull() ? hardestLightParton->pdgId() : 0 );
      else if( hadronFlavour!=0 && std::abs(partonFlavour)!=hadronFlavour )
        partonFlavour = hadronFlavour;
    }
 }

Member Data Documentation

const edm::EDGetTokenT<reco::GenParticleRefVector> JetFlavourClustering::bHadronsToken_

private

Definition at line 182 of file JetFlavourClustering.cc.

Referenced by produce().

const edm::EDGetTokenT<reco::GenParticleRefVector> JetFlavourClustering::cHadronsToken_

private

Definition at line 183 of file JetFlavourClustering.cc.

Referenced by produce().

ClusterSequencePtr JetFlavourClustering::fjClusterSeq_

private

Definition at line 196 of file JetFlavourClustering.cc.

Referenced by produce().

JetDefPtr JetFlavourClustering::fjJetDefinition_

private

Definition at line 197 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

const double JetFlavourClustering::ghostRescaling_

private

Definition at line 190 of file JetFlavourClustering.cc.

Referenced by produce().

edm::EDGetTokenT<edm::View<reco::Jet> > JetFlavourClustering::groomedJetsToken_

private

Definition at line 180 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

const bool JetFlavourClustering::hadronFlavourHasPriority_

private

Definition at line 192 of file JetFlavourClustering.cc.

Referenced by setFlavours().

const std::string JetFlavourClustering::jetAlgorithm_

private

Definition at line 187 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering().

const double JetFlavourClustering::jetPtMin_

private

Definition at line 189 of file JetFlavourClustering.cc.

Referenced by produce().

const edm::EDGetTokenT<edm::View<reco::Jet> > JetFlavourClustering::jetsToken_

private

Definition at line 179 of file JetFlavourClustering.cc.

Referenced by produce().

edm::EDGetTokenT<reco::GenParticleRefVector> JetFlavourClustering::leptonsToken_

private

Definition at line 185 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

const edm::EDGetTokenT<reco::GenParticleRefVector> JetFlavourClustering::partonsToken_

private

Definition at line 184 of file JetFlavourClustering.cc.

Referenced by produce().

const double JetFlavourClustering::relPtTolerance_

private

Definition at line 191 of file JetFlavourClustering.cc.

Referenced by produce().

const double JetFlavourClustering::rParam_

private

Definition at line 188 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), matchGroomedJets(), and matchReclusteredJets().

edm::EDGetTokenT<edm::View<reco::Jet> > JetFlavourClustering::subjetsToken_

private

Definition at line 181 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

const bool JetFlavourClustering::useLeptons_

private

Definition at line 194 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

const bool JetFlavourClustering::useSubjets_

private

Definition at line 193 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation