CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes
JetFlavourClustering Class Reference

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

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

Inheritance diagram for JetFlavourClustering:
edm::EDProducer edm::ProducerBase edm::EDConsumerBase edm::ProductRegistryHelper

Public Member Functions

 JetFlavourClustering (const edm::ParameterSet &)
 
 ~JetFlavourClustering ()
 
- Public Member Functions inherited from edm::EDProducer
 EDProducer ()
 
ModuleDescription const & moduleDescription () const
 
virtual ~EDProducer ()
 
- Public Member Functions inherited from edm::ProducerBase
 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
 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
 
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::EDProducer
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)
 
virtual void beginJob ()
 
virtual void beginLuminosityBlock (edm::LuminosityBlock &, edm::EventSetup const &)
 
virtual void beginRun (edm::Run &, edm::EventSetup const &)
 
virtual void endJob ()
 
virtual void endLuminosityBlock (edm::LuminosityBlock &, edm::EventSetup const &)
 
virtual void endRun (edm::Run &, edm::EventSetup const &)
 
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_
 
const 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_
 
const edm::EDGetTokenT
< reco::GenParticleRefVector
leptonsToken_
 
const edm::EDGetTokenT
< reco::GenParticleRefVector
partonsToken_
 
const double rParam_
 
const edm::EDGetTokenT
< edm::View< reco::Jet > > 
subjetsToken_
 
const bool useLeptons_
 
const bool useSubjets_
 

Additional Inherited Members

- Public Types inherited from edm::EDProducer
typedef EDProducer ModuleType
 
- Public Types inherited from edm::ProducerBase
typedef
ProductRegistryHelper::TypeLabelList 
TypeLabelList
 
- Protected Member Functions inherited from edm::ProducerBase
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
- 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 JetFlavorInfo 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:

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

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 which is enabled by default. The priority is given to the hadron-based flavour as follows:

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

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 flavor 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 reclustering the constituents of groomed fat jets will generally result in a jet collection different from the input groomed fat jets.

In addition, "ghost" leptons can also be clustered inside jets but they are not used in any way to determine the jet flavor. This functionality is disabled by default.

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

Definition at line 139 of file JetFlavourClustering.cc.

Constructor & Destructor Documentation

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

Definition at line 211 of file JetFlavourClustering.cc.

References edm::hlt::Exception, fjJetDefinition_, jetAlgorithm_, rParam_, and useSubjets_.

211  :
212 
214  groomedJetsToken_(mayConsume<edm::View<reco::Jet> >( iConfig.exists("groomedJets") ? iConfig.getParameter<edm::InputTag>("groomedJets") : edm::InputTag() )),
215  subjetsToken_(mayConsume<edm::View<reco::Jet> >( iConfig.exists("subjets") ? iConfig.getParameter<edm::InputTag>("subjets") : edm::InputTag() )),
216  bHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("bHadrons") )),
217  cHadronsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("cHadrons") )),
218  partonsToken_(consumes<reco::GenParticleRefVector>( iConfig.getParameter<edm::InputTag>("partons") )),
219  leptonsToken_(mayConsume<reco::GenParticleRefVector>( iConfig.exists("leptons") ? iConfig.getParameter<edm::InputTag>("leptons") : edm::InputTag() )),
220  jetAlgorithm_(iConfig.getParameter<std::string>("jetAlgorithm")),
221  rParam_(iConfig.getParameter<double>("rParam")),
222  jetPtMin_(0.), // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
223  ghostRescaling_(iConfig.exists("ghostRescaling") ? iConfig.getParameter<double>("ghostRescaling") : 1e-18),
224  hadronFlavourHasPriority_(iConfig.getParameter<bool>("hadronFlavourHasPriority")),
225  useSubjets_(iConfig.exists("groomedJets") && iConfig.exists("subjets")),
226  useLeptons_(iConfig.exists("leptons"))
227 
228 {
229  // register your products
230  produces<reco::JetFlavourInfoMatchingCollection>();
231  if( useSubjets_ )
232  produces<reco::JetFlavourInfoMatchingCollection>("SubJets");
233 
234  // set jet algorithm
235  if (jetAlgorithm_=="Kt")
236  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::kt_algorithm, rParam_) );
237  else if (jetAlgorithm_=="CambridgeAachen")
238  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam_) );
239  else if (jetAlgorithm_=="AntiKt")
240  fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam_) );
241  else
242  throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm_ << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
243 }
T getParameter(std::string const &) const
boost::shared_ptr< fastjet::JetDefinition > JetDefPtr
bool exists(std::string const &parameterName) const
checks if a parameter exists
const edm::EDGetTokenT< reco::GenParticleRefVector > bHadronsToken_
const edm::EDGetTokenT< reco::GenParticleRefVector > partonsToken_
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
const edm::EDGetTokenT< reco::GenParticleRefVector > cHadronsToken_
const std::string jetAlgorithm_
const edm::EDGetTokenT< edm::View< reco::Jet > > groomedJetsToken_
const edm::EDGetTokenT< edm::View< reco::Jet > > jetsToken_
const edm::EDGetTokenT< reco::GenParticleRefVector > leptonsToken_
EDGetTokenT< ProductType > mayConsume(edm::InputTag const &tag)
const edm::EDGetTokenT< edm::View< reco::Jet > > subjetsToken_
JetFlavourClustering::~JetFlavourClustering ( )

Definition at line 246 of file JetFlavourClustering.cc.

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

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 639 of file JetFlavourClustering.cc.

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

Referenced by produce().

643 {
644  // loop over clustered particles and assign them to different subjets based on smallest dR
645  for(reco::GenParticleRefVector::const_iterator it = clusteredParticles.begin(); it != clusteredParticles.end(); ++it)
646  {
647  std::vector<double> dR2toSubjets;
648 
649  for(size_t sj=0; sj<subjetIndices.size(); ++sj)
650  dR2toSubjets.push_back( reco::deltaR2( (*it)->rapidity(), (*it)->phi(), subjets->at(subjetIndices.at(sj)).rapidity(), subjets->at(subjetIndices.at(sj)).phi() ) );
651 
652  // find the closest subjet
653  int closestSubjetIdx = std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
654 
655  assignedParticles.at(closestSubjetIdx).push_back( *it );
656  }
657 }
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:249
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:244
const_reference at(size_type pos) const
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
Definition: deltaR.h:58
Definition: DDAxes.h:10
void JetFlavourClustering::beginJob ( void  )
privatevirtual

Reimplemented from edm::EDProducer.

Definition at line 661 of file JetFlavourClustering.cc.

662 {
663 }
void JetFlavourClustering::beginLuminosityBlock ( edm::LuminosityBlock ,
edm::EventSetup const &   
)
privatevirtual

Definition at line 684 of file JetFlavourClustering.cc.

685 {
686 }
void JetFlavourClustering::beginRun ( edm::Run ,
edm::EventSetup const &   
)
privatevirtual

Definition at line 672 of file JetFlavourClustering.cc.

673 {
674 }
void JetFlavourClustering::endJob ( void  )
privatevirtual

Reimplemented from edm::EDProducer.

Definition at line 667 of file JetFlavourClustering.cc.

667  {
668 }
void JetFlavourClustering::endLuminosityBlock ( edm::LuminosityBlock ,
edm::EventSetup const &   
)
privatevirtual

Definition at line 690 of file JetFlavourClustering.cc.

691 {
692 }
void JetFlavourClustering::endRun ( edm::Run ,
edm::EventSetup const &   
)
privatevirtual

Definition at line 678 of file JetFlavourClustering.cc.

679 {
680 }
void JetFlavourClustering::fillDescriptions ( edm::ConfigurationDescriptions descriptions)
static

Definition at line 696 of file JetFlavourClustering.cc.

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

696  {
697  //The following says we do not know what parameters are allowed so do no validation
698  // Please change this to state exactly what you do use, even if it is no parameters
700  desc.setUnknown();
701  descriptions.addDefault(desc);
702 }
void addDefault(ParameterSetDescription const &psetDescription)
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 451 of file JetFlavourClustering.cc.

References AlCaHLTBitMon_ParallelJobs::p.

Referenced by produce().

455 {
456  // insert "ghost" particles in the vector of jet constituents
457  for(reco::GenParticleRefVector::const_iterator it = particles->begin(); it != particles->end(); ++it)
458  {
459  fastjet::PseudoJet p((*it)->px(),(*it)->py(),(*it)->pz(),(*it)->energy());
460  p*=ghostRescaling; // rescale particle momentum
461  p.set_user_info(new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
462  constituents.push_back(p);
463  }
464 }
bool isLepton(const Candidate &part)
Definition: pdgIdUtils.h:19
bool isParton(const reco::Candidate &c)
Definition: CandMCTag.cc:48
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 501 of file JetFlavourClustering.cc.

References reco::deltaR(), spr::find(), j, and fwrapper::jets.

Referenced by produce().

504 {
505  std::vector<bool> jetLocks(jets->size(),false);
506  std::vector<int> jetIndices;
507 
508  for(size_t gj=0; gj<groomedJets->size(); ++gj)
509  {
510  double matchedDR = 1e9;
511  int matchedIdx = -1;
512 
513  for(size_t j=0; j<jets->size(); ++j)
514  {
515  if( jetLocks.at(j) ) continue; // skip jets that have already been matched
516 
517  double tempDR = reco::deltaR( jets->at(j).rapidity(), jets->at(j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi() );
518  if( tempDR < matchedDR )
519  {
520  matchedDR = tempDR;
521  matchedIdx = j;
522  }
523  }
524 
525  if( matchedIdx>=0 ) jetLocks.at(matchedIdx) = true;
526  jetIndices.push_back(matchedIdx);
527  }
528 
529  if( std::find( jetIndices.begin(), jetIndices.end(), -1 ) != jetIndices.end() )
530  edm::LogError("JetMatchingFailed") << "Matching groomed to original jets failed. Please check that the two jet collections belong to each other.";
531 
532  for(size_t j=0; j<jets->size(); ++j)
533  {
534  std::vector<int>::iterator matchedIndex = std::find( jetIndices.begin(), jetIndices.end(), j );
535 
536  matchedIndices.push_back( matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(),matchedIndex) : -1 );
537  }
538 }
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
int j
Definition: DBlmapReader.cc:9
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
size_type size() const
const_reference at(size_type pos) const
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 468 of file JetFlavourClustering.cc.

References reco::deltaR(), spr::find(), j, and fwrapper::jets.

Referenced by produce().

471 {
472  std::vector<bool> matchedLocks(reclusteredJets.size(),false);
473 
474  for(size_t j=0; j<jets->size(); ++j)
475  {
476  double matchedDR = 1e9;
477  int matchedIdx = -1;
478 
479  for(size_t rj=0; rj<reclusteredJets.size(); ++rj)
480  {
481  if( matchedLocks.at(rj) ) continue; // skip jets that have already been matched
482 
483  double tempDR = reco::deltaR( jets->at(j).rapidity(), jets->at(j).phi(), reclusteredJets.at(rj).rapidity(), reclusteredJets.at(rj).phi_std() );
484  if( tempDR < matchedDR )
485  {
486  matchedDR = tempDR;
487  matchedIdx = rj;
488  }
489  }
490 
491  if( matchedIdx>=0 ) matchedLocks.at(matchedIdx) = true;
492  matchedIndices.push_back(matchedIdx);
493  }
494 
495  if( std::find( matchedIndices.begin(), matchedIndices.end(), -1 ) != matchedIndices.end() )
496  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.";
497 }
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
int j
Definition: DBlmapReader.cc:9
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
size_type size() const
const_reference at(size_type pos) const
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 542 of file JetFlavourClustering.cc.

References g, and alignCSCRings::s.

Referenced by produce().

546 {
547  for(size_t g=0; g<groomedIndices.size(); ++g)
548  {
549  std::vector<int> subjetIndices;
550 
551  if( groomedIndices.at(g)>=0 )
552  {
553  for(size_t s=0; s<groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s)
554  {
555  const edm::Ptr<reco::Candidate> & subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
556 
557  for(size_t sj=0; sj<subjets->size(); ++sj)
558  {
559  if( subjet == edm::Ptr<reco::Candidate>(subjets->ptrAt(sj)) )
560  {
561  subjetIndices.push_back(sj);
562  break;
563  }
564  }
565  }
566 
567  if( subjetIndices.size() == 0 )
568  edm::LogError("SubjetMatchingFailed") << "Matching subjets to original jets failed. Please check that the groomed jet and subjet collections belong to each other.";
569 
570  matchedIndices.push_back(subjetIndices);
571  }
572  else
573  matchedIndices.push_back(subjetIndices);
574  }
575 }
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
Ptr< value_type > ptrAt(size_type i) const
size_type size() const
const_reference at(size_type pos) const
void JetFlavourClustering::produce ( edm::Event iEvent,
const edm::EventSetup iSetup 
)
privatevirtual

Implements edm::EDProducer.

Definition at line 261 of file JetFlavourClustering.cc.

References assignToSubjets(), bHadronsToken_, cHadronsToken_, alignCSCRings::e, fjClusterSeq_, fjJetDefinition_, edm::Event::getByToken(), ghostRescaling_, groomedJetsToken_, i, insertGhosts(), GhostInfo::isbHadron(), GhostInfo::isHadron(), GhostInfo::isLepton(), GhostInfo::isParton(), jetPtMin_, fwrapper::jets, jetsToken_, EgammaValidation_Wenu_cff::leptons, leptonsToken_, m, matchGroomedJets(), matchReclusteredJets(), matchSubjets(), GhostInfo::particleRef(), partonsToken_, RecoTauCleanerPlugins::pt, edm::RefVector< C, T, F >::push_back(), edm::Event::put(), setFlavours(), subjetsToken_, useLeptons_, and useSubjets_.

262 {
264  iEvent.getByToken(jetsToken_, jets);
265 
266  edm::Handle<edm::View<reco::Jet> > groomedJets;
268  if( useSubjets_ )
269  {
270  iEvent.getByToken(groomedJetsToken_, groomedJets);
271  iEvent.getByToken(subjetsToken_, subjets);
272  }
273 
275  iEvent.getByToken(bHadronsToken_, bHadrons);
276 
278  iEvent.getByToken(cHadronsToken_, cHadrons);
279 
281  iEvent.getByToken(partonsToken_, partons);
282 
284  if( useLeptons_ )
285  iEvent.getByToken(leptonsToken_, leptons);
286 
287  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> jetFlavourInfos( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(jets)) );
288  std::auto_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
289  if( useSubjets_ )
290  subjetFlavourInfos = std::auto_ptr<reco::JetFlavourInfoMatchingCollection>( new reco::JetFlavourInfoMatchingCollection(reco::JetRefBaseProd(subjets)) );
291 
292  // vector of constituents for reclustering jets and "ghosts"
293  std::vector<fastjet::PseudoJet> fjInputs;
294  // loop over all input jets and collect all their constituents
295  for(edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); ++it)
296  {
297  std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
298  std::vector<edm::Ptr<reco::Candidate> >::const_iterator m;
299  for( m = constituents.begin(); m != constituents.end(); ++m )
300  {
301  reco::CandidatePtr constit = *m;
302  if(constit->pt() == 0)
303  {
304  edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
305  continue;
306  }
307  fjInputs.push_back(fastjet::PseudoJet(constit->px(),constit->py(),constit->pz(),constit->energy()));
308  }
309  }
310  // insert "ghost" b hadrons in the vector of constituents
311  insertGhosts(bHadrons, ghostRescaling_, true, true, false, false, fjInputs);
312  // insert "ghost" c hadrons in the vector of constituents
313  insertGhosts(cHadrons, ghostRescaling_, true, false, false, false, fjInputs);
314  // insert "ghost" partons in the vector of constituents
315  insertGhosts(partons, ghostRescaling_, false, false, true, false, fjInputs);
316  // if used, insert "ghost" leptons in the vector of constituents
317  if( useLeptons_ )
318  insertGhosts(leptons, ghostRescaling_, false, false, false, true, fjInputs);
319 
320  // define jet clustering sequence
321  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs, *fjJetDefinition_ ) );
322  // recluster jet constituents and inserted "ghosts"
323  std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt( fjClusterSeq_->inclusive_jets(jetPtMin_) );
324 
325  if( inclusiveJets.size() < jets->size() )
326  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.";
327 
328  // match reclustered and original jets
329  std::vector<int> reclusteredIndices;
330  matchReclusteredJets(jets,inclusiveJets,reclusteredIndices);
331 
332  // match groomed and original jets
333  std::vector<int> groomedIndices;
334  if( useSubjets_ )
335  {
336  if( groomedJets->size() > jets->size() )
337  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.";
338 
339  matchGroomedJets(jets,groomedJets,groomedIndices);
340  }
341 
342  // match subjets and original jets
343  std::vector<std::vector<int> > subjetIndices;
344  if( useSubjets_ )
345  {
346  matchSubjets(groomedIndices,groomedJets,subjets,subjetIndices);
347  }
348 
349  // determine jet flavour
350  for(size_t i=0; i<jets->size(); ++i)
351  {
352  // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
353  if( ( fabs( inclusiveJets.at(reclusteredIndices.at(i)).pt() - jets->at(i).pt() ) / jets->at(i).pt() ) > 1e-3 ) // 0.1% difference in Pt should be sufficient to detect possible misconfigurations
354  {
355  if( jets->at(i).pt() < 10. ) // special handling for low-Pt jets (Pt<10 GeV)
356  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"
357  << "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"
358  << "Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
359  << "\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.";
360  else
361  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"
362  << "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"
363  << "\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.";
364  }
365 
366  reco::GenParticleRefVector clusteredbHadrons;
367  reco::GenParticleRefVector clusteredcHadrons;
368  reco::GenParticleRefVector clusteredPartons;
369  reco::GenParticleRefVector clusteredLeptons;
370 
371  // get jet constituents (sorted by Pt)
372  std::vector<fastjet::PseudoJet> constituents = fastjet::sorted_by_pt( inclusiveJets.at(reclusteredIndices.at(i)).constituents() );
373 
374  // loop over jet constituents and try to find "ghosts"
375  for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
376  {
377  if( !it->has_user_info() ) continue; // skip if not a "ghost"
378 
379  // "ghost" hadron
380  if( it->user_info<GhostInfo>().isHadron() )
381  {
382  // "ghost" b hadron
383  if( it->user_info<GhostInfo>().isbHadron() )
384  clusteredbHadrons.push_back(it->user_info<GhostInfo>().particleRef());
385  // "ghost" c hadron
386  else
387  clusteredcHadrons.push_back(it->user_info<GhostInfo>().particleRef());
388  }
389  // "ghost" parton
390  else if( it->user_info<GhostInfo>().isParton() )
391  clusteredPartons.push_back(it->user_info<GhostInfo>().particleRef());
392  // "ghost" lepton
393  else if( it->user_info<GhostInfo>().isLepton() )
394  clusteredLeptons.push_back(it->user_info<GhostInfo>().particleRef());
395  }
396 
397  int hadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
398  int partonFlavour = 0; // default parton flavour set to 0 (= undefined)
399 
400  // set hadron- and parton-based flavours
401  setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
402 
403  // set the JetFlavourInfo for this jet
404  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
405 
406  // if subjets are used, determine their flavour
407  if( useSubjets_ )
408  {
409  if( subjetIndices.at(i).size()==0 ) continue; // continue if the original jet does not have subjets assigned
410 
411  // define vectors of GenParticleRefVectors for hadrons and partons assigned to different subjets
412  std::vector<reco::GenParticleRefVector> assignedbHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
413  std::vector<reco::GenParticleRefVector> assignedcHadrons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
414  std::vector<reco::GenParticleRefVector> assignedPartons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
415  std::vector<reco::GenParticleRefVector> assignedLeptons(subjetIndices.at(i).size(),reco::GenParticleRefVector());
416 
417  // loop over clustered b hadrons and assign them to different subjets based on smallest dR
418  assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(i), assignedbHadrons);
419  // loop over clustered c hadrons and assign them to different subjets based on smallest dR
420  assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(i), assignedcHadrons);
421  // loop over clustered partons and assign them to different subjets based on smallest dR
422  assignToSubjets(clusteredPartons, subjets, subjetIndices.at(i), assignedPartons);
423  // if used, loop over clustered leptons and assign them to different subjets based on smallest dR
424  if( useLeptons_ )
425  assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(i), assignedLeptons);
426 
427  // loop over subjets and determine their flavour
428  for(size_t sj=0; sj<subjetIndices.at(i).size(); ++sj)
429  {
430  int subjetHadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
431  int subjetPartonFlavour = 0; // default parton flavour set to 0 (= undefined)
432 
433  // set hadron- and parton-based flavours
434  setFlavours(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
435 
436  // set the JetFlavourInfo for this subjet
437  (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] = reco::JetFlavourInfo(assignedbHadrons.at(sj), assignedcHadrons.at(sj), assignedPartons.at(sj), assignedLeptons.at(sj), subjetHadronFlavour, subjetPartonFlavour);
438  }
439  }
440  }
441 
442  // put jet flavour infos in the event
443  iEvent.put( jetFlavourInfos );
444  // put subjet flavour infos in the event
445  if( useSubjets_ )
446  iEvent.put( subjetFlavourInfos, "SubJets" );
447 }
const reco::GenParticleRef & particleRef() const
int i
Definition: DBlmapReader.cc:9
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
void assignToSubjets(const reco::GenParticleRefVector &clusteredParticles, const edm::Handle< edm::View< reco::Jet > > &subjets, const std::vector< int > &subjetIndices, std::vector< reco::GenParticleRefVector > &assignedParticles)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
boost::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
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 setFlavours(const reco::GenParticleRefVector &clusteredbHadrons, const reco::GenParticleRefVector &clusteredcHadrons, const reco::GenParticleRefVector &clusteredPartons, int &hadronFlavour, int &partonFlavour)
const edm::EDGetTokenT< reco::GenParticleRefVector > bHadronsToken_
void matchGroomedJets(const edm::Handle< edm::View< reco::Jet > > &jets, const edm::Handle< edm::View< reco::Jet > > &matchedJets, std::vector< int > &matchedIndices)
const edm::EDGetTokenT< reco::GenParticleRefVector > partonsToken_
const bool isParton() const
const bool isLepton() const
const edm::EDGetTokenT< reco::GenParticleRefVector > cHadronsToken_
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:116
vector< PseudoJet > jets
Class storing the jet flavour information.
void matchReclusteredJets(const edm::Handle< edm::View< reco::Jet > > &jets, const std::vector< fastjet::PseudoJet > &matchedJets, std::vector< int > &matchedIndices)
const edm::EDGetTokenT< edm::View< reco::Jet > > groomedJetsToken_
edm::RefVector< GenParticleCollection > GenParticleRefVector
vector of reference to GenParticle in the same collection
const edm::EDGetTokenT< edm::View< reco::Jet > > jetsToken_
const bool isbHadron() const
const bool isHadron() const
void push_back(value_type const &ref)
Add a Ref&lt;C, T&gt; to the RefVector.
Definition: RefVector.h:64
const edm::EDGetTokenT< reco::GenParticleRefVector > leptonsToken_
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)
ClusterSequencePtr fjClusterSeq_
const edm::EDGetTokenT< edm::View< reco::Jet > > subjetsToken_
void JetFlavourClustering::setFlavours ( const reco::GenParticleRefVector clusteredbHadrons,
const reco::GenParticleRefVector clusteredcHadrons,
const reco::GenParticleRefVector clusteredPartons,
int &  hadronFlavour,
int &  partonFlavour 
)
private

Definition at line 579 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().

584 {
585  reco::GenParticleRef hardestParton;
586  reco::GenParticleRef hardestLightParton;
587  reco::GenParticleRef flavourParton;
588 
589  // loop over clustered partons (already sorted by Pt)
590  for(reco::GenParticleRefVector::const_iterator it = clusteredPartons.begin(); it != clusteredPartons.end(); ++it)
591  {
592  // hardest parton
593  if( hardestParton.isNull() )
594  hardestParton = (*it);
595  // hardest light-flavour parton
596  if( hardestLightParton.isNull() )
597  {
598  if( CandMCTagUtils::isLightParton( *(*it) ) )
599  hardestLightParton = (*it);
600  }
601  // c flavour
602  if( flavourParton.isNull() && ( abs( (*it)->pdgId() ) == 4 ) )
603  flavourParton = (*it);
604  // b flavour gets priority
605  if( abs( (*it)->pdgId() ) == 5 )
606  {
607  if( flavourParton.isNull() )
608  flavourParton = (*it);
609  else if( abs( flavourParton->pdgId() ) != 5 )
610  flavourParton = (*it);
611  }
612  }
613 
614  // set hadron-based flavour
615  if( clusteredbHadrons.size()>0 )
616  hadronFlavour = 5;
617  else if( clusteredcHadrons.size()>0 && clusteredbHadrons.size()==0 )
618  hadronFlavour = 4;
619  // set parton-based flavour
620  if( flavourParton.isNull() )
621  {
622  if( hardestParton.isNonnull() ) partonFlavour = hardestParton->pdgId();
623  }
624  else
625  partonFlavour = flavourParton->pdgId();
626 
627  // if enabled, check for conflicts between hadron- and parton-based flavours and give priority to the hadron-based flavour
629  {
630  if( hadronFlavour==0 && (abs(partonFlavour)==4 || abs(partonFlavour)==5) )
631  partonFlavour = ( hardestLightParton.isNonnull() ? hardestLightParton->pdgId() : 0 );
632  else if( hadronFlavour!=0 && abs(partonFlavour)!=hadronFlavour )
633  partonFlavour = hadronFlavour;
634  }
635 }
bool isLightParton(const reco::Candidate &c)
Definition: CandMCTag.cc:63
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:249
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:244
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
bool isNull() const
Checks for null.
Definition: Ref.h:247
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
size_type size() const
Size of the RefVector.
Definition: RefVector.h:89

Member Data Documentation

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

Definition at line 187 of file JetFlavourClustering.cc.

Referenced by produce().

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

Definition at line 188 of file JetFlavourClustering.cc.

Referenced by produce().

ClusterSequencePtr JetFlavourClustering::fjClusterSeq_
private

Definition at line 200 of file JetFlavourClustering.cc.

Referenced by produce().

JetDefPtr JetFlavourClustering::fjJetDefinition_
private

Definition at line 201 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

const double JetFlavourClustering::ghostRescaling_
private

Definition at line 195 of file JetFlavourClustering.cc.

Referenced by produce().

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

Definition at line 185 of file JetFlavourClustering.cc.

Referenced by produce().

const bool JetFlavourClustering::hadronFlavourHasPriority_
private

Definition at line 196 of file JetFlavourClustering.cc.

Referenced by setFlavours().

const std::string JetFlavourClustering::jetAlgorithm_
private

Definition at line 192 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering().

const double JetFlavourClustering::jetPtMin_
private

Definition at line 194 of file JetFlavourClustering.cc.

Referenced by produce().

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

Definition at line 184 of file JetFlavourClustering.cc.

Referenced by produce().

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

Definition at line 190 of file JetFlavourClustering.cc.

Referenced by produce().

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

Definition at line 189 of file JetFlavourClustering.cc.

Referenced by produce().

const double JetFlavourClustering::rParam_
private

Definition at line 193 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering().

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

Definition at line 186 of file JetFlavourClustering.cc.

Referenced by produce().

const bool JetFlavourClustering::useLeptons_
private

Definition at line 198 of file JetFlavourClustering.cc.

Referenced by produce().

const bool JetFlavourClustering::useSubjets_
private

Definition at line 197 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().