CMS 3D CMS Logo

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::stream::EDProducer<>

Public Member Functions

 JetFlavourClustering (const edm::ParameterSet &)
 
 ~JetFlavourClustering () override
 
- Public Member Functions inherited from edm::stream::EDProducer<>
 EDProducer ()=default
 
 EDProducer (const EDProducer &)=delete
 
bool hasAbilityToProduceInBeginLumis () const final
 
bool hasAbilityToProduceInBeginProcessBlocks () const final
 
bool hasAbilityToProduceInBeginRuns () const final
 
bool hasAbilityToProduceInEndLumis () const final
 
bool hasAbilityToProduceInEndProcessBlocks () const final
 
bool hasAbilityToProduceInEndRuns () const final
 
const EDProduceroperator= (const EDProducer &)=delete
 

Static Public Member Functions

static void fillDescriptions (edm::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)
 
void produce (edm::Event &, const edm::EventSetup &) override
 
void setFlavours (const reco::GenParticleRefVector &clusteredbHadrons, const reco::GenParticleRefVector &clusteredcHadrons, const reco::GenParticleRefVector &clusteredPartons, int &hadronFlavour, int &partonFlavour)
 

Private Attributes

const edm::EDGetTokenT< reco::GenParticleRefVectorbHadronsToken_
 
const edm::EDGetTokenT< reco::GenParticleRefVectorcHadronsToken_
 
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::GenParticleRefVectorleptonsToken_
 
const edm::EDGetTokenT< reco::GenParticleRefVectorpartonsToken_
 
const double relPtTolerance_
 
const double rParam_
 
edm::EDGetTokenT< edm::View< reco::Jet > > subjetsToken_
 
const bool useLeptons_
 
const bool useSubjets_
 
edm::EDGetTokenT< edm::ValueMap< float > > weightsToken_
 

Additional Inherited Members

- Public Types inherited from edm::stream::EDProducer<>
using CacheTypes = CacheContexts< T... >
 
using GlobalCache = typename CacheTypes::GlobalCache
 
using HasAbility = AbilityChecker< T... >
 
using InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache
 
using LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache
 
using LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >
 
using LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache
 
using RunCache = typename CacheTypes::RunCache
 
using RunContext = RunContextT< RunCache, GlobalCache >
 
using RunSummaryCache = typename CacheTypes::RunSummaryCache
 

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:

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

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:

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

Constructor & Destructor Documentation

◆ JetFlavourClustering()

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

Definition at line 215 of file JetFlavourClustering.cc.

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

216  : jetsToken_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
217  bHadronsToken_(consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("bHadrons"))),
218  cHadronsToken_(consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("cHadrons"))),
219  partonsToken_(consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("partons"))),
220  jetAlgorithm_(iConfig.getParameter<std::string>("jetAlgorithm")),
221  rParam_(iConfig.getParameter<double>("rParam")),
222 
223  jetPtMin_(
224  0.), // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
225  ghostRescaling_(iConfig.exists("ghostRescaling") ? iConfig.getParameter<double>("ghostRescaling") : 1e-18),
227  iConfig.exists("relPtTolerance")
228  ? iConfig.getParameter<double>("relPtTolerance")
229  : 1e-03), // 0.1% relative difference in Pt should be sufficient to detect possible misconfigurations
230  hadronFlavourHasPriority_(iConfig.getParameter<bool>("hadronFlavourHasPriority")),
231  useSubjets_(iConfig.exists("groomedJets") && iConfig.exists("subjets")),
232  useLeptons_(iConfig.exists("leptons"))
233 
234 {
235  // register your products
236  produces<reco::JetFlavourInfoMatchingCollection>();
237  if (iConfig.existsAs<edm::InputTag>("weights"))
238  weightsToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("weights"));
239 
240  if (useSubjets_)
241  produces<reco::JetFlavourInfoMatchingCollection>("SubJets");
242 
243  // set jet algorithm
244  if (jetAlgorithm_ == "Kt")
245  fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::kt_algorithm, rParam_);
246  else if (jetAlgorithm_ == "CambridgeAachen")
247  fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::cambridge_algorithm, rParam_);
248  else if (jetAlgorithm_ == "AntiKt")
249  fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
250  else
251  throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm_
252  << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
253 
254  if (useSubjets_) {
255  groomedJetsToken_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("groomedJets"));
256  subjetsToken_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("subjets"));
257  }
258  if (useLeptons_) {
259  leptonsToken_ = consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("leptons"));
260  }
261 }
edm::EDGetTokenT< reco::GenParticleRefVector > leptonsToken_
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
edm::EDGetTokenT< edm::View< reco::Jet > > subjetsToken_
bool exists(std::string const &parameterName) const
checks if a parameter exists
const edm::EDGetTokenT< reco::GenParticleRefVector > bHadronsToken_
bool existsAs(std::string const &parameterName, bool trackiness=true) const
checks if a parameter exists as a given type
Definition: ParameterSet.h:172
edm::EDGetTokenT< edm::ValueMap< float > > weightsToken_
const edm::EDGetTokenT< reco::GenParticleRefVector > partonsToken_
const edm::EDGetTokenT< reco::GenParticleRefVector > cHadronsToken_
const std::string jetAlgorithm_
const edm::EDGetTokenT< edm::View< reco::Jet > > jetsToken_
edm::EDGetTokenT< edm::View< reco::Jet > > groomedJetsToken_

◆ ~JetFlavourClustering()

JetFlavourClustering::~JetFlavourClustering ( )
override

Definition at line 263 of file JetFlavourClustering.cc.

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

Member Function Documentation

◆ assignToSubjets()

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

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

Referenced by produce().

734  {
735  // loop over clustered particles and assign them to different subjets based on smallest dR
736  for (reco::GenParticleRefVector::const_iterator it = clusteredParticles.begin(); it != clusteredParticles.end();
737  ++it) {
738  std::vector<double> dR2toSubjets;
739 
740  for (size_t sj = 0; sj < subjetIndices.size(); ++sj)
741  dR2toSubjets.push_back(reco::deltaR2((*it)->rapidity(),
742  (*it)->phi(),
743  subjets->at(subjetIndices.at(sj)).rapidity(),
744  subjets->at(subjetIndices.at(sj)).phi()));
745 
746  // find the closest subjet
747  int closestSubjetIdx =
748  std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
749 
750  assignedParticles.at(closestSubjetIdx).push_back(*it);
751  }
752 }
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:228
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:223

◆ fillDescriptions()

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

Definition at line 755 of file JetFlavourClustering.cc.

References edm::ConfigurationDescriptions::addDefault(), and submitPVResolutionJobs::desc.

755  {
756  //The following says we do not know what parameters are allowed so do no validation
757  // Please change this to state exactly what you do use, even if it is no parameters
759  desc.setUnknown();
760  descriptions.addDefault(desc);
761 }
void addDefault(ParameterSetDescription const &psetDescription)

◆ insertGhosts()

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

References AK4PFJetsMCFlavourInfos_cfi::ghostRescaling, reco::isLepton(), CandMCTagUtils::isParton(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, AlCaHLTBitMon_ParallelJobs::p, and ecalTrigSettings_cff::particles.

Referenced by produce().

544  {
545  // insert "ghost" particles in the vector of jet constituents
546  for (reco::GenParticleRefVector::const_iterator it = particles->begin(); it != particles->end(); ++it) {
547  if ((*it)->pt() == 0) {
548  edm::LogInfo("NullTransverseMomentum") << "dropping input ghost candidate with pt=0";
549  continue;
550  }
551  fastjet::PseudoJet p((*it)->px(), (*it)->py(), (*it)->pz(), (*it)->energy());
552  p *= ghostRescaling; // rescale particle momentum
553  p.set_user_info(new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
554  constituents.push_back(p);
555  }
556 }
bool isLepton(const Candidate &part)
Definition: pdgIdUtils.h:13
Log< level::Info, false > LogInfo
bool isParton(const reco::Candidate &c)
Definition: CandMCTag.cc:46

◆ matchGroomedJets()

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

References reco::deltaR2(), HLT_2024v13_cff::distance, spr::find(), dqmiolumiharvest::j, PDWG_EXODelayedJetMET_cff::jets, rParam_, and mathSSE::sqrt().

Referenced by produce().

602  {
603  std::vector<bool> jetLocks(jets->size(), false);
604  std::vector<int> jetIndices;
605 
606  for (size_t gj = 0; gj < groomedJets->size(); ++gj) {
607  double matchedDR2 = 1e9;
608  int matchedIdx = -1;
609 
610  if (groomedJets->at(gj).pt() > 0.) // skip pathological cases of groomed jets with Pt=0
611  {
612  for (size_t j = 0; j < jets->size(); ++j) {
613  if (jetLocks.at(j))
614  continue; // skip jets that have already been matched
615 
616  double tempDR2 = reco::deltaR2(
617  jets->at(j).rapidity(), jets->at(j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi());
618  if (tempDR2 < matchedDR2) {
619  matchedDR2 = tempDR2;
620  matchedIdx = j;
621  }
622  }
623  }
624 
625  if (matchedIdx >= 0) {
626  if (matchedDR2 > rParam_ * rParam_) {
627  edm::LogWarning("MatchedJetsFarApart")
628  << "Matched groomed jet " << gj << " and original jet " << matchedIdx
629  << " are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam_
630  << ".\n"
631  << "This is not expected so the matching of these two jets has been discarded. Please check that the two "
632  "jet collections belong to each other.";
633  matchedIdx = -1;
634  } else
635  jetLocks.at(matchedIdx) = true;
636  }
637  jetIndices.push_back(matchedIdx);
638  }
639 
640  for (size_t j = 0; j < jets->size(); ++j) {
641  std::vector<int>::iterator matchedIndex = std::find(jetIndices.begin(), jetIndices.end(), j);
642 
643  matchedIndices.push_back(matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(), matchedIndex) : -1);
644  }
645 }
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
T sqrt(T t)
Definition: SSEVec.h:19
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
Log< level::Warning, false > LogWarning

◆ matchReclusteredJets()

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

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

Referenced by produce().

561  {
562  std::vector<bool> matchedLocks(reclusteredJets.size(), false);
563 
564  for (size_t j = 0; j < jets->size(); ++j) {
565  double matchedDR2 = 1e9;
566  int matchedIdx = -1;
567 
568  for (size_t rj = 0; rj < reclusteredJets.size(); ++rj) {
569  if (matchedLocks.at(rj))
570  continue; // skip jets that have already been matched
571 
572  double tempDR2 = reco::deltaR2(jets->at(j).rapidity(),
573  jets->at(j).phi(),
574  reclusteredJets.at(rj).rapidity(),
575  reclusteredJets.at(rj).phi_std());
576  if (tempDR2 < matchedDR2) {
577  matchedDR2 = tempDR2;
578  matchedIdx = rj;
579  }
580  }
581 
582  if (matchedIdx >= 0) {
583  if (matchedDR2 > rParam_ * rParam_) {
584  edm::LogError("JetMatchingFailed") << "Matched reclustered jet " << matchedIdx << " and original jet " << j
585  << " are separated by dR=" << sqrt(matchedDR2)
586  << " which is greater than the jet size R=" << rParam_ << ".\n"
587  << "This is not expected so please check that the jet algorithm and jet "
588  "size match those used for the original jet collection.";
589  } else
590  matchedLocks.at(matchedIdx) = true;
591  } else
592  edm::LogError("JetMatchingFailed") << "Matching reclustered to original jets failed. Please check that the jet "
593  "algorithm and jet size match those used for the original jet collection.";
594 
595  matchedIndices.push_back(matchedIdx);
596  }
597 }
Log< level::Error, false > LogError
T sqrt(T t)
Definition: SSEVec.h:19
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16

◆ matchSubjets()

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

References g, and alignCSCRings::s.

Referenced by produce().

651  {
652  for (size_t g = 0; g < groomedIndices.size(); ++g) {
653  std::vector<int> subjetIndices;
654 
655  if (groomedIndices.at(g) >= 0) {
656  for (size_t s = 0; s < groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s) {
657  const edm::Ptr<reco::Candidate>& subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
658 
659  for (size_t sj = 0; sj < subjets->size(); ++sj) {
660  if (subjet == edm::Ptr<reco::Candidate>(subjets->ptrAt(sj))) {
661  subjetIndices.push_back(sj);
662  break;
663  }
664  }
665  }
666 
667  if (subjetIndices.empty())
668  edm::LogError("SubjetMatchingFailed") << "Matching subjets to original jets failed. Please check that the "
669  "groomed jet and subjet collections belong to each other.";
670 
671  matchedIndices.push_back(subjetIndices);
672  } else
673  matchedIndices.push_back(subjetIndices);
674  }
675 }
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

◆ produce()

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

Definition at line 273 of file JetFlavourClustering.cc.

References funct::abs(), assignToSubjets(), AK4GenJetFlavourInfos_cfi::bHadrons, bHadronsToken_, AK4GenJetFlavourInfos_cfi::cHadrons, cHadronsToken_, fjClusterSeq_, fjJetDefinition_, ghostRescaling_, groomedJetsToken_, custom_jme_cff::hadronFlavour, mps_fire::i, iEvent, insertGhosts(), edm::Ptr< T >::isAvailable(), GhostInfo::isbHadron(), GhostInfo::isHadron(), GhostInfo::isLepton(), edm::Ptr< T >::isNonnull(), GhostInfo::isParton(), edm::EDGetTokenT< T >::isUninitialized(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, GenHFHadronMatcher_cfi::jetFlavourInfos, jetPtMin_, PDWG_EXODelayedJetMET_cff::jets, jetsToken_, HLT_2024v13_cff::leptons, leptonsToken_, visualization-live-secondInstance_cfg::m, matchGroomedJets(), matchReclusteredJets(), matchSubjets(), eostools::move(), GhostInfo::particleRef(), custom_jme_cff::partonFlavour, dqmAnalyzer_cff::partons, partonsToken_, DiDispStaMuonMonitor_cfi::pt, edm::RefVector< C, T, F >::push_back(), relPtTolerance_, setFlavours(), subjetsToken_, useLeptons_, useSubjets_, w(), hltDeepSecondaryVertexTagInfosPFPuppi_cfi::weights, and weightsToken_.

273  {
275  iEvent.getByToken(jetsToken_, jets);
276 
279  if (useSubjets_) {
280  iEvent.getByToken(groomedJetsToken_, groomedJets);
281  iEvent.getByToken(subjetsToken_, subjets);
282  }
283 
285  iEvent.getByToken(bHadronsToken_, bHadrons);
286 
288  iEvent.getByToken(cHadronsToken_, cHadrons);
289 
291  iEvent.getByToken(partonsToken_, partons);
292 
295  iEvent.getByToken(weightsToken_, weights);
296 
298  if (useLeptons_)
299  iEvent.getByToken(leptonsToken_, leptons);
300 
301  auto jetFlavourInfos = std::make_unique<reco::JetFlavourInfoMatchingCollection>(reco::JetRefBaseProd(jets));
302  std::unique_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
303  if (useSubjets_)
304  subjetFlavourInfos = std::make_unique<reco::JetFlavourInfoMatchingCollection>(reco::JetRefBaseProd(subjets));
305 
306  // vector of constituents for reclustering jets and "ghosts"
307  std::vector<fastjet::PseudoJet> fjInputs;
308  unsigned int reserve = jets->size() * 128 + bHadrons->size() + cHadrons->size() + partons->size();
309  if (useLeptons_)
310  reserve += leptons->size();
311  fjInputs.reserve(reserve);
312  // loop over all input jets and collect all their constituents
313  for (edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); ++it) {
314  std::vector<edm::Ptr<reco::Candidate>> constituents = it->getJetConstituents();
315  std::vector<edm::Ptr<reco::Candidate>>::const_iterator m;
316  for (m = constituents.begin(); m != constituents.end(); ++m) {
317  const reco::CandidatePtr& constit = *m;
318  if (!constit.isNonnull() || !constit.isAvailable()) {
319  edm::LogError("MissingJetConstituent") << "Jet constituent required for jet reclustering is missing. "
320  "Reclustered jets are not guaranteed to reproduce the original jets!";
321  continue;
322  }
323  if (constit->pt() == 0) {
324  edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
325  continue;
326  }
327  if (it->isWeighted()) {
329  throw cms::Exception("MissingConstituentWeight")
330  << "JetFlavourClustering: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
331  float w = (*weights)[constit];
332  fjInputs.push_back(
333  fastjet::PseudoJet(constit->px() * w, constit->py() * w, constit->pz() * w, constit->energy() * w));
334  } else {
335  fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
336  }
337  }
338  }
339  // insert "ghost" b hadrons in the vector of constituents
340  insertGhosts(bHadrons, ghostRescaling_, true, true, false, false, fjInputs);
341  // insert "ghost" c hadrons in the vector of constituents
342  insertGhosts(cHadrons, ghostRescaling_, true, false, false, false, fjInputs);
343  // insert "ghost" partons in the vector of constituents
344  insertGhosts(partons, ghostRescaling_, false, false, true, false, fjInputs);
345  // if used, insert "ghost" leptons in the vector of constituents
346  if (useLeptons_)
347  insertGhosts(leptons, ghostRescaling_, false, false, false, true, fjInputs);
348 
349  // define jet clustering sequence
350  fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs, *fjJetDefinition_);
351  // recluster jet constituents and inserted "ghosts"
352  std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
353 
354  if (inclusiveJets.size() < jets->size())
355  edm::LogError("TooFewReclusteredJets")
356  << "There are fewer reclustered (" << inclusiveJets.size() << ") than original jets (" << jets->size()
357  << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
358 
359  // match reclustered and original jets
360  std::vector<int> reclusteredIndices;
361  matchReclusteredJets(jets, inclusiveJets, reclusteredIndices);
362 
363  // match groomed and original jets
364  std::vector<int> groomedIndices;
365  if (useSubjets_) {
366  if (groomedJets->size() > jets->size())
367  edm::LogError("TooManyGroomedJets")
368  << "There are more groomed (" << groomedJets->size() << ") than original jets (" << jets->size()
369  << "). Please check that the two jet collections belong to each other.";
370 
371  matchGroomedJets(jets, groomedJets, groomedIndices);
372  }
373 
374  // match subjets and original jets
375  std::vector<std::vector<int>> subjetIndices;
376  if (useSubjets_) {
377  matchSubjets(groomedIndices, groomedJets, subjets, subjetIndices);
378  }
379 
380  // determine jet flavour
381  for (size_t i = 0; i < jets->size(); ++i) {
382  reco::GenParticleRefVector clusteredbHadrons;
383  reco::GenParticleRefVector clusteredcHadrons;
384  reco::GenParticleRefVector clusteredPartons;
385  reco::GenParticleRefVector clusteredLeptons;
386 
387  // if matching reclustered to original jets failed
388  if (reclusteredIndices.at(i) < 0) {
389  // set an empty JetFlavourInfo for this jet
390  (*jetFlavourInfos)[jets->refAt(i)] =
391  reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
392  } else if (jets->at(i).pt() == 0) {
393  edm::LogWarning("NullTransverseMomentum")
394  << "The original jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
395 
396  // set an empty JetFlavourInfo for this jet
397  (*jetFlavourInfos)[jets->refAt(i)] =
398  reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
399 
400  // if subjets are used
401  if (useSubjets_ && !subjetIndices.at(i).empty()) {
402  // loop over subjets
403  for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj) {
404  // set an empty JetFlavourInfo for this subjet
405  (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] =
410  0,
411  0);
412  }
413  }
414  } else {
415  // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
416  if ((std::abs(inclusiveJets.at(reclusteredIndices.at(i)).pt() - jets->at(i).pt()) / jets->at(i).pt()) >
417  relPtTolerance_) {
418  if (jets->at(i).pt() < 10.) // special handling for low-Pt jets (Pt<10 GeV)
419  edm::LogWarning("JetPtMismatchAtLowPt")
420  << "The reclustered and original jet " << i << " have different Pt's ("
421  << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << jets->at(i).pt()
422  << " GeV, respectively).\n"
423  << "Please check that the jet algorithm and jet size match those used for the original jet collection "
424  "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
425  "CaloJets which are presently not supported.\n"
426  << "Since the mismatch is at low Pt (Pt<10 GeV), it is ignored and only a warning is issued.\n"
427  << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision "
428  "in which case make sure the original jet collection is produced and reclustering is performed in the "
429  "same job.";
430  else
431  edm::LogError("JetPtMismatch")
432  << "The reclustered and original jet " << i << " have different Pt's ("
433  << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << jets->at(i).pt()
434  << " GeV, respectively).\n"
435  << "Please check that the jet algorithm and jet size match those used for the original jet collection "
436  "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
437  "CaloJets which are presently not supported.\n"
438  << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision "
439  "in which case make sure the original jet collection is produced and reclustering is performed in the "
440  "same job.";
441  }
442 
443  // get jet constituents (sorted by Pt)
444  std::vector<fastjet::PseudoJet> constituents =
445  fastjet::sorted_by_pt(inclusiveJets.at(reclusteredIndices.at(i)).constituents());
446 
447  // loop over jet constituents and try to find "ghosts"
448  for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it) {
449  if (!it->has_user_info())
450  continue; // skip if not a "ghost"
451 
452  // "ghost" hadron
453  if (it->user_info<GhostInfo>().isHadron()) {
454  // "ghost" b hadron
455  if (it->user_info<GhostInfo>().isbHadron())
456  clusteredbHadrons.push_back(it->user_info<GhostInfo>().particleRef());
457  // "ghost" c hadron
458  else
459  clusteredcHadrons.push_back(it->user_info<GhostInfo>().particleRef());
460  }
461  // "ghost" parton
462  else if (it->user_info<GhostInfo>().isParton())
463  clusteredPartons.push_back(it->user_info<GhostInfo>().particleRef());
464  // "ghost" lepton
465  else if (it->user_info<GhostInfo>().isLepton())
466  clusteredLeptons.push_back(it->user_info<GhostInfo>().particleRef());
467  }
468 
469  int hadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
470  int partonFlavour = 0; // default parton flavour set to 0 (= undefined)
471 
472  // set hadron- and parton-based flavours
473  setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
474 
475  // set the JetFlavourInfo for this jet
476  (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(
477  clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
478  }
479 
480  // if subjets are used, determine their flavour
481  if (useSubjets_) {
482  if (subjetIndices.at(i).empty())
483  continue; // continue if the original jet does not have subjets assigned
484 
485  // define vectors of GenParticleRefVectors for hadrons and partons assigned to different subjets
486  std::vector<reco::GenParticleRefVector> assignedbHadrons(subjetIndices.at(i).size(),
488  std::vector<reco::GenParticleRefVector> assignedcHadrons(subjetIndices.at(i).size(),
490  std::vector<reco::GenParticleRefVector> assignedPartons(subjetIndices.at(i).size(), reco::GenParticleRefVector());
491  std::vector<reco::GenParticleRefVector> assignedLeptons(subjetIndices.at(i).size(), reco::GenParticleRefVector());
492 
493  // loop over clustered b hadrons and assign them to different subjets based on smallest dR
494  assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(i), assignedbHadrons);
495  // loop over clustered c hadrons and assign them to different subjets based on smallest dR
496  assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(i), assignedcHadrons);
497  // loop over clustered partons and assign them to different subjets based on smallest dR
498  assignToSubjets(clusteredPartons, subjets, subjetIndices.at(i), assignedPartons);
499  // if used, loop over clustered leptons and assign them to different subjets based on smallest dR
500  if (useLeptons_)
501  assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(i), assignedLeptons);
502 
503  // loop over subjets and determine their flavour
504  for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj) {
505  int subjetHadronFlavour = 0; // default hadron flavour set to 0 (= undefined)
506  int subjetPartonFlavour = 0; // default parton flavour set to 0 (= undefined)
507 
508  // set hadron- and parton-based flavours
509  setFlavours(assignedbHadrons.at(sj),
510  assignedcHadrons.at(sj),
511  assignedPartons.at(sj),
512  subjetHadronFlavour,
513  subjetPartonFlavour);
514 
515  // set the JetFlavourInfo for this subjet
516  (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] =
517  reco::JetFlavourInfo(assignedbHadrons.at(sj),
518  assignedcHadrons.at(sj),
519  assignedPartons.at(sj),
520  assignedLeptons.at(sj),
521  subjetHadronFlavour,
522  subjetPartonFlavour);
523  }
524  }
525  }
526 
527  //deallocate only at the end of the event processing
528  fjClusterSeq_.reset();
529 
530  // put jet flavour infos in the event
532  // put subjet flavour infos in the event
533  if (useSubjets_)
534  iEvent.put(std::move(subjetFlavourInfos), "SubJets");
535 }
edm::EDGetTokenT< reco::GenParticleRefVector > leptonsToken_
const bool isLepton() const
T w() 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)
const bool isParton() const
edm::EDGetTokenT< edm::View< reco::Jet > > subjetsToken_
void assignToSubjets(const reco::GenParticleRefVector &clusteredParticles, const edm::Handle< edm::View< reco::Jet >> &subjets, const std::vector< int > &subjetIndices, std::vector< reco::GenParticleRefVector > &assignedParticles)
constexpr bool isUninitialized() const noexcept
Definition: EDGetToken.h:104
void setFlavours(const reco::GenParticleRefVector &clusteredbHadrons, const reco::GenParticleRefVector &clusteredcHadrons, const reco::GenParticleRefVector &clusteredPartons, int &hadronFlavour, int &partonFlavour)
const edm::EDGetTokenT< reco::GenParticleRefVector > bHadronsToken_
Log< level::Error, false > LogError
edm::EDGetTokenT< edm::ValueMap< float > > weightsToken_
bool isAvailable() const
Definition: Ptr.h:232
const edm::EDGetTokenT< reco::GenParticleRefVector > partonsToken_
edm::RefToBaseProd< reco::Jet > JetRefBaseProd
Definition: JetCollection.h:13
const edm::EDGetTokenT< reco::GenParticleRefVector > cHadronsToken_
int iEvent
Definition: GenABIO.cc:224
void matchReclusteredJets(const edm::Handle< edm::View< reco::Jet >> &jets, const std::vector< fastjet::PseudoJet > &matchedJets, std::vector< int > &matchedIndices)
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
Class storing the jet flavour information.
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void matchGroomedJets(const edm::Handle< edm::View< reco::Jet >> &jets, const edm::Handle< edm::View< reco::Jet >> &matchedJets, std::vector< int > &matchedIndices)
edm::RefVector< GenParticleCollection > GenParticleRefVector
vector of reference to GenParticle in the same collection
const reco::GenParticleRef & particleRef() const
const edm::EDGetTokenT< edm::View< reco::Jet > > jetsToken_
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:88
edm::EDGetTokenT< edm::View< reco::Jet > > groomedJetsToken_
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67
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)
Log< level::Warning, false > LogWarning
const bool isbHadron() const
ClusterSequencePtr fjClusterSeq_
def move(src, dest)
Definition: eostools.py:511
const bool isHadron() const

◆ setFlavours()

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

Definition at line 678 of file JetFlavourClustering.cc.

References funct::abs(), edm::RefVector< C, T, F >::begin(), edm::RefVector< C, T, F >::empty(), edm::RefVector< C, T, F >::end(), custom_jme_cff::hadronFlavour, hadronFlavourHasPriority_, CandMCTagUtils::isLightParton(), edm::Ref< C, T, F >::isNonnull(), edm::Ref< C, T, F >::isNull(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, and custom_jme_cff::partonFlavour.

Referenced by produce().

682  {
683  reco::GenParticleRef hardestParton;
684  reco::GenParticleRef hardestLightParton;
685  reco::GenParticleRef flavourParton;
686 
687  // loop over clustered partons (already sorted by Pt)
688  for (reco::GenParticleRefVector::const_iterator it = clusteredPartons.begin(); it != clusteredPartons.end(); ++it) {
689  // hardest parton
690  if (hardestParton.isNull())
691  hardestParton = (*it);
692  // hardest light-flavour parton
693  if (hardestLightParton.isNull()) {
695  hardestLightParton = (*it);
696  }
697  // c flavour
698  if (flavourParton.isNull() && (std::abs((*it)->pdgId()) == 4))
699  flavourParton = (*it);
700  // b flavour gets priority
701  if (std::abs((*it)->pdgId()) == 5) {
702  if (flavourParton.isNull())
703  flavourParton = (*it);
704  else if (std::abs(flavourParton->pdgId()) != 5)
705  flavourParton = (*it);
706  }
707  }
708 
709  // set hadron-based flavour
710  if (!clusteredbHadrons.empty())
711  hadronFlavour = 5;
712  else if (!clusteredcHadrons.empty() && clusteredbHadrons.empty())
713  hadronFlavour = 4;
714  // set parton-based flavour
715  if (flavourParton.isNull()) {
716  if (hardestParton.isNonnull())
717  partonFlavour = hardestParton->pdgId();
718  } else
719  partonFlavour = flavourParton->pdgId();
720 
721  // if enabled, check for conflicts between hadron- and parton-based flavours and give priority to the hadron-based flavour
723  if (hadronFlavour == 0 && (std::abs(partonFlavour) == 4 || std::abs(partonFlavour) == 5))
724  partonFlavour = (hardestLightParton.isNonnull() ? hardestLightParton->pdgId() : 0);
725  else if (hadronFlavour != 0 && std::abs(partonFlavour) != hadronFlavour)
727  }
728 }
bool isLightParton(const reco::Candidate &c)
Definition: CandMCTag.cc:55
bool empty() const
Is the RefVector empty.
Definition: RefVector.h:99
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool isNull() const
Checks for null.
Definition: Ref.h:235
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:228
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:223

Member Data Documentation

◆ bHadronsToken_

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

Definition at line 188 of file JetFlavourClustering.cc.

Referenced by produce().

◆ cHadronsToken_

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

Definition at line 189 of file JetFlavourClustering.cc.

Referenced by produce().

◆ fjClusterSeq_

ClusterSequencePtr JetFlavourClustering::fjClusterSeq_
private

Definition at line 204 of file JetFlavourClustering.cc.

Referenced by produce().

◆ fjJetDefinition_

JetDefPtr JetFlavourClustering::fjJetDefinition_
private

Definition at line 205 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

◆ ghostRescaling_

const double JetFlavourClustering::ghostRescaling_
private

Definition at line 197 of file JetFlavourClustering.cc.

Referenced by produce().

◆ groomedJetsToken_

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

Definition at line 186 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

◆ hadronFlavourHasPriority_

const bool JetFlavourClustering::hadronFlavourHasPriority_
private

Definition at line 199 of file JetFlavourClustering.cc.

Referenced by setFlavours().

◆ jetAlgorithm_

const std::string JetFlavourClustering::jetAlgorithm_
private

Definition at line 194 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering().

◆ jetPtMin_

const double JetFlavourClustering::jetPtMin_
private

Definition at line 196 of file JetFlavourClustering.cc.

Referenced by produce().

◆ jetsToken_

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

Definition at line 185 of file JetFlavourClustering.cc.

Referenced by produce().

◆ leptonsToken_

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

Definition at line 192 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

◆ partonsToken_

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

Definition at line 190 of file JetFlavourClustering.cc.

Referenced by produce().

◆ relPtTolerance_

const double JetFlavourClustering::relPtTolerance_
private

Definition at line 198 of file JetFlavourClustering.cc.

Referenced by produce().

◆ rParam_

const double JetFlavourClustering::rParam_
private

◆ subjetsToken_

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

Definition at line 187 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

◆ useLeptons_

const bool JetFlavourClustering::useLeptons_
private

Definition at line 202 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

◆ useSubjets_

const bool JetFlavourClustering::useSubjets_
private

Definition at line 200 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().

◆ weightsToken_

edm::EDGetTokenT<edm::ValueMap<float> > JetFlavourClustering::weightsToken_
private

Definition at line 191 of file JetFlavourClustering.cc.

Referenced by JetFlavourClustering(), and produce().