JetMCAlgos/plugins/JetFlavourClustering.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    JetFlavourClustering
0004 // Class:      JetFlavourClustering
0005 //
0006 /**\class JetFlavourClustering JetFlavourClustering.cc PhysicsTools/JetMCAlgos/plugins/JetFlavourClustering.cc
0007  * \brief Clusters hadrons, partons, and jet contituents to determine the jet flavour
0008  *
0009  * This producer clusters hadrons, partons and jet contituents to determine the jet flavour. The jet flavour information
0010  * is stored in the event as an AssociationVector which associates an object of type JetFlavourInfo to each of the jets.
0011  *
0012  * The producer takes as input jets and hadron and partons selected by the HadronAndPartonSelector producer. The hadron
0013  * and parton four-momenta are rescaled by a very small number (default rescale factor is 10e-18) which turns them into
0014  * the so-called "ghosts". The "ghost" hadrons and partons are clustered together with all of the jet constituents. It is
0015  * important to use the same clustering algorithm and jet size as for the original input jet collection. Since the
0016  * "ghost" hadrons and partons are extremely soft, the resulting jet collection will be practically identical to the
0017  * original one but now with "ghost" hadrons and partons clustered inside jets. The jet flavour is determined based on
0018  * the "ghost" hadrons clustered inside a jet:
0019  *
0020  * - jet is considered a b jet if there is at least one b "ghost" hadron clustered inside it (hadronFlavour = 5)
0021  *
0022  * - jet is considered a c jet if there is at least one c and no b "ghost" hadrons clustered inside it (hadronFlavour = 4)
0023  *
0024  * - jet is considered a light-flavour jet if there are no b or c "ghost" hadrons clustered inside it (hadronFlavour = 0)
0025  *
0026  * To further assign a more specific flavour to light-flavour jets, "ghost" partons are used:
0027  *
0028  * - jet is considered a b jet if there is at least one b "ghost" parton clustered inside it (partonFlavour = 5)
0029  *
0030  * - jet is considered a c jet if there is at least one c and no b "ghost" partons clustered inside it (partonFlavour = 4)
0031  *
0032  * - jet is considered a light-flavour jet if there are light-flavour and no b or c "ghost" partons clustered inside it.
0033  *   The jet is assigned the flavour of the hardest light-flavour "ghost" parton clustered inside it (partonFlavour = 1, 2, 3, or 21)
0034  *
0035  * - jet has an undefined flavour if there are no "ghost" partons clustered inside it (partonFlavour = 0)
0036  *
0037  * In rare instances a conflict between the hadron- and parton-based flavours can occur. In such cases it is possible to
0038  * keep both flavours or to give priority to the hadron-based flavour. This is controlled by the 'hadronFlavourHasPriority'
0039  * switch. The priority is given to the hadron-based flavour as follows:
0040  *
0041  * - if hadronFlavour==0 && (partonFlavour==4 || partonFlavour==5): partonFlavour is set to the flavour of the hardest
0042  *   light-flavour parton clustered inside the jet if such parton exists. Otherwise, the parton flavour is left undefined
0043  *
0044  * - if hadronFlavour!=0 && hadronFlavour!=partonFlavour: partonFlavour is set equal to hadronFlavour
0045  *
0046  * The producer is also capable of assigning the flavour to subjets of fat jets, in which case it produces an additional
0047  * AssociationVector providing the flavour information for subjets. In order to assign the flavour to subjets, three input
0048  * jet collections are required:
0049  *
0050  * - jets, in this case represented by fat jets
0051  *
0052  * - 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)
0053  *
0054  * - subjets, derived from the groomed fat jets
0055  *
0056  * The "ghost" hadrons and partons clustered inside a fat jet are assigned to the closest subjet in the rapidity-phi
0057  * space. Once hadrons and partons have been assigned to subjets, the subjet flavour is determined in the same way as for
0058  * jets. The reason for requiring three jet collections as input in order to determine the subjet flavour is to avoid
0059  * possible inconsistencies between the fat jet and subjet flavours (such as a non-b fat jet having a b subjet and vice
0060  * versa) as well as the fact that re-clustering the constituents of groomed fat jets will generally result in a jet
0061  * collection different from the input groomed fat jets. Also note that "ghost" particles generally cannot be clustered
0062  * inside subjets in the same way this is done for fat jets. This is because some of the jet grooming techniques could
0063  * reject such very soft particle. So instead, the "ghost" particles are assigned to the closest subjet.
0064  *
0065  * Finally, "ghost" leptons can also be clustered inside jets but they are not used in any way to determine the jet
0066  * flavour. This functionality is optional and is potentially useful to identify jets from hadronic taus.
0067  *
0068  * For more details, please refer to
0069  * https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideBTagMCTools
0070  *
0071  */
0072 //
0073 // Original Author:  Dinko Ferencek
0074 //         Created:  Wed Nov  6 00:49:55 CET 2013
0075 //
0076 //
0077
0078 // system include files
0079 #include <memory>
0080
0081 // user include files
0082 #include "FWCore/Framework/interface/Frameworkfwd.h"
0083 #include "FWCore/Framework/interface/stream/EDProducer.h"
0084
0085 #include "FWCore/Framework/interface/Event.h"
0086 #include "FWCore/Framework/interface/MakerMacros.h"
0087
0088 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0089
0090 #include "DataFormats/JetReco/interface/Jet.h"
0091 #include "DataFormats/JetReco/interface/JetCollection.h"
0092 #include "DataFormats/JetMatching/interface/JetFlavourInfo.h"
0093 #include "DataFormats/JetMatching/interface/JetFlavourInfoMatching.h"
0094 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0095 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0096 #include "DataFormats/Math/interface/deltaR.h"
0097 #include "PhysicsTools/JetMCUtils/interface/CandMCTag.h"
0098 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0099
0100 #include "fastjet/JetDefinition.hh"
0101 #include "fastjet/ClusterSequence.hh"
0102 #include "fastjet/Selector.hh"
0103 #include "fastjet/PseudoJet.hh"
0104
0105 //
0106 // constants, enums and typedefs
0107 //
0108 typedef std::shared_ptr<fastjet::ClusterSequence> ClusterSequencePtr;
0109 typedef std::shared_ptr<fastjet::JetDefinition> JetDefPtr;
0110
0111 //
0112 // class declaration
0113 //
0114 class GhostInfo : public fastjet::PseudoJet::UserInfoBase {
0115 public:
0116   GhostInfo(const bool& isHadron,
0117             const bool& isbHadron,
0118             const bool& isParton,
0119             const bool& isLepton,
0120             const reco::GenParticleRef& particleRef)
0121       : m_particleRef(particleRef) {
0122     m_type = 0;
0123     if (isHadron)
0124       m_type |= (1 << 0);
0125     if (isbHadron)
0126       m_type |= (1 << 1);
0127     if (isParton)
0128       m_type |= (1 << 2);
0129     if (isLepton)
0130       m_type |= (1 << 3);
0131   }
0132
0133   const bool isHadron() const { return (m_type & (1 << 0)); }
0134   const bool isbHadron() const { return (m_type & (1 << 1)); }
0135   const bool isParton() const { return (m_type & (1 << 2)); }
0136   const bool isLepton() const { return (m_type & (1 << 3)); }
0137   const reco::GenParticleRef& particleRef() const { return m_particleRef; }
0138
0139 protected:
0140   const reco::GenParticleRef m_particleRef;
0141   int m_type;
0142 };
0143
0144 class JetFlavourClustering : public edm::stream::EDProducer<> {
0145 public:
0146   explicit JetFlavourClustering(const edm::ParameterSet&);
0147   ~JetFlavourClustering() override;
0148
0149   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0150
0151 private:
0152   void produce(edm::Event&, const edm::EventSetup&) override;
0153
0154   void insertGhosts(const edm::Handle<reco::GenParticleRefVector>& particles,
0155                     const double ghostRescaling,
0156                     const bool isHadron,
0157                     const bool isbHadron,
0158                     const bool isParton,
0159                     const bool isLepton,
0160                     std::vector<fastjet::PseudoJet>& constituents);
0161
0162   void matchReclusteredJets(const edm::Handle<edm::View<reco::Jet>>& jets,
0163                             const std::vector<fastjet::PseudoJet>& matchedJets,
0164                             std::vector<int>& matchedIndices);
0165   void matchGroomedJets(const edm::Handle<edm::View<reco::Jet>>& jets,
0166                         const edm::Handle<edm::View<reco::Jet>>& matchedJets,
0167                         std::vector<int>& matchedIndices);
0168   void matchSubjets(const std::vector<int>& groomedIndices,
0169                     const edm::Handle<edm::View<reco::Jet>>& groomedJets,
0170                     const edm::Handle<edm::View<reco::Jet>>& subjets,
0171                     std::vector<std::vector<int>>& matchedIndices);
0172
0173   void setFlavours(const reco::GenParticleRefVector& clusteredbHadrons,
0174                    const reco::GenParticleRefVector& clusteredcHadrons,
0175                    const reco::GenParticleRefVector& clusteredPartons,
0176                    int& hadronFlavour,
0177                    int& partonFlavour);
0178
0179   void assignToSubjets(const reco::GenParticleRefVector& clusteredParticles,
0180                        const edm::Handle<edm::View<reco::Jet>>& subjets,
0181                        const std::vector<int>& subjetIndices,
0182                        std::vector<reco::GenParticleRefVector>& assignedParticles);
0183
0184   // ----------member data ---------------------------
0185   const edm::EDGetTokenT<edm::View<reco::Jet>> jetsToken_;            // Input jet collection
0186   edm::EDGetTokenT<edm::View<reco::Jet>> groomedJetsToken_;           // Input groomed jet collection
0187   edm::EDGetTokenT<edm::View<reco::Jet>> subjetsToken_;               // Input subjet collection
0188   const edm::EDGetTokenT<reco::GenParticleRefVector> bHadronsToken_;  // Input b hadron collection
0189   const edm::EDGetTokenT<reco::GenParticleRefVector> cHadronsToken_;  // Input c hadron collection
0190   const edm::EDGetTokenT<reco::GenParticleRefVector> partonsToken_;   // Input parton collection
0191   edm::EDGetTokenT<edm::ValueMap<float>> weightsToken_;               // Input weights collection
0192   edm::EDGetTokenT<reco::GenParticleRefVector> leptonsToken_;         // Input lepton collection
0193
0194   const std::string jetAlgorithm_;
0195   const double rParam_;
0196   const double jetPtMin_;
0197   const double ghostRescaling_;
0198   const double relPtTolerance_;
0199   const bool hadronFlavourHasPriority_;
0200   const bool useSubjets_;
0201
0202   const bool useLeptons_;
0203
0204   ClusterSequencePtr fjClusterSeq_;
0205   JetDefPtr fjJetDefinition_;
0206 };
0207
0208 //
0209 // static data member definitions
0210 //
0211
0212 //
0213 // constructors and destructor
0214 //
0215 JetFlavourClustering::JetFlavourClustering(const edm::ParameterSet& iConfig)
0216     : jetsToken_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
0217       bHadronsToken_(consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("bHadrons"))),
0218       cHadronsToken_(consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("cHadrons"))),
0219       partonsToken_(consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("partons"))),
0220       jetAlgorithm_(iConfig.getParameter<std::string>("jetAlgorithm")),
0221       rParam_(iConfig.getParameter<double>("rParam")),
0222
0223       jetPtMin_(
0224           0.),  // hardcoded to 0. since we simply want to recluster all input jets which already had some PtMin applied
0225       ghostRescaling_(iConfig.exists("ghostRescaling") ? iConfig.getParameter<double>("ghostRescaling") : 1e-18),
0226       relPtTolerance_(
0227           iConfig.exists("relPtTolerance")
0228               ? iConfig.getParameter<double>("relPtTolerance")
0229               : 1e-03),  // 0.1% relative difference in Pt should be sufficient to detect possible misconfigurations
0230       hadronFlavourHasPriority_(iConfig.getParameter<bool>("hadronFlavourHasPriority")),
0231       useSubjets_(iConfig.exists("groomedJets") && iConfig.exists("subjets")),
0232       useLeptons_(iConfig.exists("leptons"))
0233
0234 {
0235   // register your products
0236   produces<reco::JetFlavourInfoMatchingCollection>();
0237   if (iConfig.existsAs<edm::InputTag>("weights"))
0238     weightsToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("weights"));
0239
0240   if (useSubjets_)
0241     produces<reco::JetFlavourInfoMatchingCollection>("SubJets");
0242
0243   // set jet algorithm
0244   if (jetAlgorithm_ == "Kt")
0245     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::kt_algorithm, rParam_);
0246   else if (jetAlgorithm_ == "CambridgeAachen")
0247     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::cambridge_algorithm, rParam_);
0248   else if (jetAlgorithm_ == "AntiKt")
0249     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
0250   else
0251     throw cms::Exception("InvalidJetAlgorithm") << "Jet clustering algorithm is invalid: " << jetAlgorithm_
0252                                                 << ", use CambridgeAachen | Kt | AntiKt" << std::endl;
0253
0254   if (useSubjets_) {
0255     groomedJetsToken_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("groomedJets"));
0256     subjetsToken_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("subjets"));
0257   }
0258   if (useLeptons_) {
0259     leptonsToken_ = consumes<reco::GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("leptons"));
0260   }
0261 }
0262
0263 JetFlavourClustering::~JetFlavourClustering() {
0264   // do anything here that needs to be done at desctruction time
0265   // (e.g. close files, deallocate resources etc.)
0266 }
0267
0268 //
0269 // member functions
0270 //
0271
0272 // ------------ method called to produce the data  ------------
0273 void JetFlavourClustering::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0274   edm::Handle<edm::View<reco::Jet>> jets;
0275   iEvent.getByToken(jetsToken_, jets);
0276
0277   edm::Handle<edm::View<reco::Jet>> groomedJets;
0278   edm::Handle<edm::View<reco::Jet>> subjets;
0279   if (useSubjets_) {
0280     iEvent.getByToken(groomedJetsToken_, groomedJets);
0281     iEvent.getByToken(subjetsToken_, subjets);
0282   }
0283
0284   edm::Handle<reco::GenParticleRefVector> bHadrons;
0285   iEvent.getByToken(bHadronsToken_, bHadrons);
0286
0287   edm::Handle<reco::GenParticleRefVector> cHadrons;
0288   iEvent.getByToken(cHadronsToken_, cHadrons);
0289
0290   edm::Handle<reco::GenParticleRefVector> partons;
0291   iEvent.getByToken(partonsToken_, partons);
0292
0293   edm::Handle<edm::ValueMap<float>> weights;
0294   if (!weightsToken_.isUninitialized())
0295     iEvent.getByToken(weightsToken_, weights);
0296
0297   edm::Handle<reco::GenParticleRefVector> leptons;
0298   if (useLeptons_)
0299     iEvent.getByToken(leptonsToken_, leptons);
0300
0301   auto jetFlavourInfos = std::make_unique<reco::JetFlavourInfoMatchingCollection>(reco::JetRefBaseProd(jets));
0302   std::unique_ptr<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfos;
0303   if (useSubjets_)
0304     subjetFlavourInfos = std::make_unique<reco::JetFlavourInfoMatchingCollection>(reco::JetRefBaseProd(subjets));
0305
0306   // vector of constituents for reclustering jets and "ghosts"
0307   std::vector<fastjet::PseudoJet> fjInputs;
0308   unsigned int reserve = jets->size() * 128 + bHadrons->size() + cHadrons->size() + partons->size();
0309   if (useLeptons_)
0310     reserve += leptons->size();
0311   fjInputs.reserve(reserve);
0312   // loop over all input jets and collect all their constituents
0313   for (edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); ++it) {
0314     std::vector<edm::Ptr<reco::Candidate>> constituents = it->getJetConstituents();
0315     std::vector<edm::Ptr<reco::Candidate>>::const_iterator m;
0316     for (m = constituents.begin(); m != constituents.end(); ++m) {
0317       const reco::CandidatePtr& constit = *m;
0318       if (!constit.isNonnull() || !constit.isAvailable()) {
0319         edm::LogError("MissingJetConstituent") << "Jet constituent required for jet reclustering is missing. "
0320                                                   "Reclustered jets are not guaranteed to reproduce the original jets!";
0321         continue;
0322       }
0323       if (constit->pt() == 0) {
0324         edm::LogWarning("NullTransverseMomentum") << "dropping input candidate with pt=0";
0325         continue;
0326       }
0327       if (it->isWeighted()) {
0328         if (weightsToken_.isUninitialized())
0329           throw cms::Exception("MissingConstituentWeight")
0330               << "JetFlavourClustering: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0331         float w = (*weights)[constit];
0332         fjInputs.push_back(
0333             fastjet::PseudoJet(constit->px() * w, constit->py() * w, constit->pz() * w, constit->energy() * w));
0334       } else {
0335         fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
0336       }
0337     }
0338   }
0339   // insert "ghost" b hadrons in the vector of constituents
0340   insertGhosts(bHadrons, ghostRescaling_, true, true, false, false, fjInputs);
0341   // insert "ghost" c hadrons in the vector of constituents
0342   insertGhosts(cHadrons, ghostRescaling_, true, false, false, false, fjInputs);
0343   // insert "ghost" partons in the vector of constituents
0344   insertGhosts(partons, ghostRescaling_, false, false, true, false, fjInputs);
0345   // if used, insert "ghost" leptons in the vector of constituents
0346   if (useLeptons_)
0347     insertGhosts(leptons, ghostRescaling_, false, false, false, true, fjInputs);
0348
0349   // define jet clustering sequence
0350   fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs, *fjJetDefinition_);
0351   // recluster jet constituents and inserted "ghosts"
0352   std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
0353
0354   if (inclusiveJets.size() < jets->size())
0355     edm::LogError("TooFewReclusteredJets")
0356         << "There are fewer reclustered (" << inclusiveJets.size() << ") than original jets (" << jets->size()
0357         << "). Please check that the jet algorithm and jet size match those used for the original jet collection.";
0358
0359   // match reclustered and original jets
0360   std::vector<int> reclusteredIndices;
0361   matchReclusteredJets(jets, inclusiveJets, reclusteredIndices);
0362
0363   // match groomed and original jets
0364   std::vector<int> groomedIndices;
0365   if (useSubjets_) {
0366     if (groomedJets->size() > jets->size())
0367       edm::LogError("TooManyGroomedJets")
0368           << "There are more groomed (" << groomedJets->size() << ") than original jets (" << jets->size()
0369           << "). Please check that the two jet collections belong to each other.";
0370
0371     matchGroomedJets(jets, groomedJets, groomedIndices);
0372   }
0373
0374   // match subjets and original jets
0375   std::vector<std::vector<int>> subjetIndices;
0376   if (useSubjets_) {
0377     matchSubjets(groomedIndices, groomedJets, subjets, subjetIndices);
0378   }
0379
0380   // determine jet flavour
0381   for (size_t i = 0; i < jets->size(); ++i) {
0382     reco::GenParticleRefVector clusteredbHadrons;
0383     reco::GenParticleRefVector clusteredcHadrons;
0384     reco::GenParticleRefVector clusteredPartons;
0385     reco::GenParticleRefVector clusteredLeptons;
0386
0387     // if matching reclustered to original jets failed
0388     if (reclusteredIndices.at(i) < 0) {
0389       // set an empty JetFlavourInfo for this jet
0390       (*jetFlavourInfos)[jets->refAt(i)] =
0391           reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
0392     } else if (jets->at(i).pt() == 0) {
0393       edm::LogWarning("NullTransverseMomentum")
0394           << "The original jet " << i << " has Pt=0. This is not expected so the jet will be skipped.";
0395
0396       // set an empty JetFlavourInfo for this jet
0397       (*jetFlavourInfos)[jets->refAt(i)] =
0398           reco::JetFlavourInfo(clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, 0, 0);
0399
0400       // if subjets are used
0401       if (useSubjets_ && !subjetIndices.at(i).empty()) {
0402         // loop over subjets
0403         for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj) {
0404           // set an empty JetFlavourInfo for this subjet
0405           (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] =
0406               reco::JetFlavourInfo(reco::GenParticleRefVector(),
0407                                    reco::GenParticleRefVector(),
0408                                    reco::GenParticleRefVector(),
0409                                    reco::GenParticleRefVector(),
0410                                    0,
0411                                    0);
0412         }
0413       }
0414     } else {
0415       // since the "ghosts" are extremely soft, the configuration and ordering of the reclustered and original jets should in principle stay the same
0416       if ((std::abs(inclusiveJets.at(reclusteredIndices.at(i)).pt() - jets->at(i).pt()) / jets->at(i).pt()) >
0417           relPtTolerance_) {
0418         if (jets->at(i).pt() < 10.)  // special handling for low-Pt jets (Pt<10 GeV)
0419           edm::LogWarning("JetPtMismatchAtLowPt")
0420               << "The reclustered and original jet " << i << " have different Pt's ("
0421               << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << jets->at(i).pt()
0422               << " GeV, respectively).\n"
0423               << "Please check that the jet algorithm and jet size match those used for the original jet collection "
0424                  "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
0425                  "CaloJets which are presently not supported.\n"
0426               << "Since the mismatch is at low Pt (Pt<10 GeV), it is ignored and only a warning is issued.\n"
0427               << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision "
0428                  "in which case make sure the original jet collection is produced and reclustering is performed in the "
0429                  "same job.";
0430         else
0431           edm::LogError("JetPtMismatch")
0432               << "The reclustered and original jet " << i << " have different Pt's ("
0433               << inclusiveJets.at(reclusteredIndices.at(i)).pt() << " vs " << jets->at(i).pt()
0434               << " GeV, respectively).\n"
0435               << "Please check that the jet algorithm and jet size match those used for the original jet collection "
0436                  "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
0437                  "CaloJets which are presently not supported.\n"
0438               << "\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision "
0439                  "in which case make sure the original jet collection is produced and reclustering is performed in the "
0440                  "same job.";
0441       }
0442
0443       // get jet constituents (sorted by Pt)
0444       std::vector<fastjet::PseudoJet> constituents =
0445           fastjet::sorted_by_pt(inclusiveJets.at(reclusteredIndices.at(i)).constituents());
0446
0447       // loop over jet constituents and try to find "ghosts"
0448       for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it) {
0449         if (!it->has_user_info())
0450           continue;  // skip if not a "ghost"
0451
0452         // "ghost" hadron
0453         if (it->user_info<GhostInfo>().isHadron()) {
0454           // "ghost" b hadron
0455           if (it->user_info<GhostInfo>().isbHadron())
0456             clusteredbHadrons.push_back(it->user_info<GhostInfo>().particleRef());
0457           // "ghost" c hadron
0458           else
0459             clusteredcHadrons.push_back(it->user_info<GhostInfo>().particleRef());
0460         }
0461         // "ghost" parton
0462         else if (it->user_info<GhostInfo>().isParton())
0463           clusteredPartons.push_back(it->user_info<GhostInfo>().particleRef());
0464         // "ghost" lepton
0465         else if (it->user_info<GhostInfo>().isLepton())
0466           clusteredLeptons.push_back(it->user_info<GhostInfo>().particleRef());
0467       }
0468
0469       int hadronFlavour = 0;  // default hadron flavour set to 0 (= undefined)
0470       int partonFlavour = 0;  // default parton flavour set to 0 (= undefined)
0471
0472       // set hadron- and parton-based flavours
0473       setFlavours(clusteredbHadrons, clusteredcHadrons, clusteredPartons, hadronFlavour, partonFlavour);
0474
0475       // set the JetFlavourInfo for this jet
0476       (*jetFlavourInfos)[jets->refAt(i)] = reco::JetFlavourInfo(
0477           clusteredbHadrons, clusteredcHadrons, clusteredPartons, clusteredLeptons, hadronFlavour, partonFlavour);
0478     }
0479
0480     // if subjets are used, determine their flavour
0481     if (useSubjets_) {
0482       if (subjetIndices.at(i).empty())
0483         continue;  // continue if the original jet does not have subjets assigned
0484
0485       // define vectors of GenParticleRefVectors for hadrons and partons assigned to different subjets
0486       std::vector<reco::GenParticleRefVector> assignedbHadrons(subjetIndices.at(i).size(),
0487                                                                reco::GenParticleRefVector());
0488       std::vector<reco::GenParticleRefVector> assignedcHadrons(subjetIndices.at(i).size(),
0489                                                                reco::GenParticleRefVector());
0490       std::vector<reco::GenParticleRefVector> assignedPartons(subjetIndices.at(i).size(), reco::GenParticleRefVector());
0491       std::vector<reco::GenParticleRefVector> assignedLeptons(subjetIndices.at(i).size(), reco::GenParticleRefVector());
0492
0493       // loop over clustered b hadrons and assign them to different subjets based on smallest dR
0494       assignToSubjets(clusteredbHadrons, subjets, subjetIndices.at(i), assignedbHadrons);
0495       // loop over clustered c hadrons and assign them to different subjets based on smallest dR
0496       assignToSubjets(clusteredcHadrons, subjets, subjetIndices.at(i), assignedcHadrons);
0497       // loop over clustered partons and assign them to different subjets based on smallest dR
0498       assignToSubjets(clusteredPartons, subjets, subjetIndices.at(i), assignedPartons);
0499       // if used, loop over clustered leptons and assign them to different subjets based on smallest dR
0500       if (useLeptons_)
0501         assignToSubjets(clusteredLeptons, subjets, subjetIndices.at(i), assignedLeptons);
0502
0503       // loop over subjets and determine their flavour
0504       for (size_t sj = 0; sj < subjetIndices.at(i).size(); ++sj) {
0505         int subjetHadronFlavour = 0;  // default hadron flavour set to 0 (= undefined)
0506         int subjetPartonFlavour = 0;  // default parton flavour set to 0 (= undefined)
0507
0508         // set hadron- and parton-based flavours
0509         setFlavours(assignedbHadrons.at(sj),
0510                     assignedcHadrons.at(sj),
0511                     assignedPartons.at(sj),
0512                     subjetHadronFlavour,
0513                     subjetPartonFlavour);
0514
0515         // set the JetFlavourInfo for this subjet
0516         (*subjetFlavourInfos)[subjets->refAt(subjetIndices.at(i).at(sj))] =
0517             reco::JetFlavourInfo(assignedbHadrons.at(sj),
0518                                  assignedcHadrons.at(sj),
0519                                  assignedPartons.at(sj),
0520                                  assignedLeptons.at(sj),
0521                                  subjetHadronFlavour,
0522                                  subjetPartonFlavour);
0523       }
0524     }
0525   }
0526
0527   //deallocate only at the end of the event processing
0528   fjClusterSeq_.reset();
0529
0530   // put jet flavour infos in the event
0531   iEvent.put(std::move(jetFlavourInfos));
0532   // put subjet flavour infos in the event
0533   if (useSubjets_)
0534     iEvent.put(std::move(subjetFlavourInfos), "SubJets");
0535 }
0536
0537 // ------------ method that inserts "ghost" particles in the vector of jet constituents ------------
0538 void JetFlavourClustering::insertGhosts(const edm::Handle<reco::GenParticleRefVector>& particles,
0539                                         const double ghostRescaling,
0540                                         const bool isHadron,
0541                                         const bool isbHadron,
0542                                         const bool isParton,
0543                                         const bool isLepton,
0544                                         std::vector<fastjet::PseudoJet>& constituents) {
0545   // insert "ghost" particles in the vector of jet constituents
0546   for (reco::GenParticleRefVector::const_iterator it = particles->begin(); it != particles->end(); ++it) {
0547     if ((*it)->pt() == 0) {
0548       edm::LogInfo("NullTransverseMomentum") << "dropping input ghost candidate with pt=0";
0549       continue;
0550     }
0551     fastjet::PseudoJet p((*it)->px(), (*it)->py(), (*it)->pz(), (*it)->energy());
0552     p *= ghostRescaling;  // rescale particle momentum
0553     p.set_user_info(new GhostInfo(isHadron, isbHadron, isParton, isLepton, *it));
0554     constituents.push_back(p);
0555   }
0556 }
0557
0558 // ------------ method that matches reclustered and original jets based on minimum dR ------------
0559 void JetFlavourClustering::matchReclusteredJets(const edm::Handle<edm::View<reco::Jet>>& jets,
0560                                                 const std::vector<fastjet::PseudoJet>& reclusteredJets,
0561                                                 std::vector<int>& matchedIndices) {
0562   std::vector<bool> matchedLocks(reclusteredJets.size(), false);
0563
0564   for (size_t j = 0; j < jets->size(); ++j) {
0565     double matchedDR2 = 1e9;
0566     int matchedIdx = -1;
0567
0568     for (size_t rj = 0; rj < reclusteredJets.size(); ++rj) {
0569       if (matchedLocks.at(rj))
0570         continue;  // skip jets that have already been matched
0571
0572       double tempDR2 = reco::deltaR2(jets->at(j).rapidity(),
0573                                      jets->at(j).phi(),
0574                                      reclusteredJets.at(rj).rapidity(),
0575                                      reclusteredJets.at(rj).phi_std());
0576       if (tempDR2 < matchedDR2) {
0577         matchedDR2 = tempDR2;
0578         matchedIdx = rj;
0579       }
0580     }
0581
0582     if (matchedIdx >= 0) {
0583       if (matchedDR2 > rParam_ * rParam_) {
0584         edm::LogError("JetMatchingFailed") << "Matched reclustered jet " << matchedIdx << " and original jet " << j
0585                                            << " are separated by dR=" << sqrt(matchedDR2)
0586                                            << " which is greater than the jet size R=" << rParam_ << ".\n"
0587                                            << "This is not expected so please check that the jet algorithm and jet "
0588                                               "size match those used for the original jet collection.";
0589       } else
0590         matchedLocks.at(matchedIdx) = true;
0591     } else
0592       edm::LogError("JetMatchingFailed") << "Matching reclustered to original jets failed. Please check that the jet "
0593                                             "algorithm and jet size match those used for the original jet collection.";
0594
0595     matchedIndices.push_back(matchedIdx);
0596   }
0597 }
0598
0599 // ------------ method that matches groomed and original jets based on minimum dR ------------
0600 void JetFlavourClustering::matchGroomedJets(const edm::Handle<edm::View<reco::Jet>>& jets,
0601                                             const edm::Handle<edm::View<reco::Jet>>& groomedJets,
0602                                             std::vector<int>& matchedIndices) {
0603   std::vector<bool> jetLocks(jets->size(), false);
0604   std::vector<int> jetIndices;
0605
0606   for (size_t gj = 0; gj < groomedJets->size(); ++gj) {
0607     double matchedDR2 = 1e9;
0608     int matchedIdx = -1;
0609
0610     if (groomedJets->at(gj).pt() > 0.)  // skip pathological cases of groomed jets with Pt=0
0611     {
0612       for (size_t j = 0; j < jets->size(); ++j) {
0613         if (jetLocks.at(j))
0614           continue;  // skip jets that have already been matched
0615
0616         double tempDR2 = reco::deltaR2(
0617             jets->at(j).rapidity(), jets->at(j).phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi());
0618         if (tempDR2 < matchedDR2) {
0619           matchedDR2 = tempDR2;
0620           matchedIdx = j;
0621         }
0622       }
0623     }
0624
0625     if (matchedIdx >= 0) {
0626       if (matchedDR2 > rParam_ * rParam_) {
0627         edm::LogWarning("MatchedJetsFarApart")
0628             << "Matched groomed jet " << gj << " and original jet " << matchedIdx
0629             << " are separated by dR=" << sqrt(matchedDR2) << " which is greater than the jet size R=" << rParam_
0630             << ".\n"
0631             << "This is not expected so the matching of these two jets has been discarded. Please check that the two "
0632                "jet collections belong to each other.";
0633         matchedIdx = -1;
0634       } else
0635         jetLocks.at(matchedIdx) = true;
0636     }
0637     jetIndices.push_back(matchedIdx);
0638   }
0639
0640   for (size_t j = 0; j < jets->size(); ++j) {
0641     std::vector<int>::iterator matchedIndex = std::find(jetIndices.begin(), jetIndices.end(), j);
0642
0643     matchedIndices.push_back(matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(), matchedIndex) : -1);
0644   }
0645 }
0646
0647 // ------------ method that matches subjets and original jets ------------
0648 void JetFlavourClustering::matchSubjets(const std::vector<int>& groomedIndices,
0649                                         const edm::Handle<edm::View<reco::Jet>>& groomedJets,
0650                                         const edm::Handle<edm::View<reco::Jet>>& subjets,
0651                                         std::vector<std::vector<int>>& matchedIndices) {
0652   for (size_t g = 0; g < groomedIndices.size(); ++g) {
0653     std::vector<int> subjetIndices;
0654
0655     if (groomedIndices.at(g) >= 0) {
0656       for (size_t s = 0; s < groomedJets->at(groomedIndices.at(g)).numberOfDaughters(); ++s) {
0657         const edm::Ptr<reco::Candidate>& subjet = groomedJets->at(groomedIndices.at(g)).daughterPtr(s);
0658
0659         for (size_t sj = 0; sj < subjets->size(); ++sj) {
0660           if (subjet == edm::Ptr<reco::Candidate>(subjets->ptrAt(sj))) {
0661             subjetIndices.push_back(sj);
0662             break;
0663           }
0664         }
0665       }
0666
0667       if (subjetIndices.empty())
0668         edm::LogError("SubjetMatchingFailed") << "Matching subjets to original jets failed. Please check that the "
0669                                                  "groomed jet and subjet collections belong to each other.";
0670
0671       matchedIndices.push_back(subjetIndices);
0672     } else
0673       matchedIndices.push_back(subjetIndices);
0674   }
0675 }
0676
0677 // ------------ method that sets hadron- and parton-based flavours ------------
0678 void JetFlavourClustering::setFlavours(const reco::GenParticleRefVector& clusteredbHadrons,
0679                                        const reco::GenParticleRefVector& clusteredcHadrons,
0680                                        const reco::GenParticleRefVector& clusteredPartons,
0681                                        int& hadronFlavour,
0682                                        int& partonFlavour) {
0683   reco::GenParticleRef hardestParton;
0684   reco::GenParticleRef hardestLightParton;
0685   reco::GenParticleRef flavourParton;
0686
0687   // loop over clustered partons (already sorted by Pt)
0688   for (reco::GenParticleRefVector::const_iterator it = clusteredPartons.begin(); it != clusteredPartons.end(); ++it) {
0689     // hardest parton
0690     if (hardestParton.isNull())
0691       hardestParton = (*it);
0692     // hardest light-flavour parton
0693     if (hardestLightParton.isNull()) {
0694       if (CandMCTagUtils::isLightParton(*(*it)))
0695         hardestLightParton = (*it);
0696     }
0697     // c flavour
0698     if (flavourParton.isNull() && (std::abs((*it)->pdgId()) == 4))
0699       flavourParton = (*it);
0700     // b flavour gets priority
0701     if (std::abs((*it)->pdgId()) == 5) {
0702       if (flavourParton.isNull())
0703         flavourParton = (*it);
0704       else if (std::abs(flavourParton->pdgId()) != 5)
0705         flavourParton = (*it);
0706     }
0707   }
0708
0709   // set hadron-based flavour
0710   if (!clusteredbHadrons.empty())
0711     hadronFlavour = 5;
0712   else if (!clusteredcHadrons.empty() && clusteredbHadrons.empty())
0713     hadronFlavour = 4;
0714   // set parton-based flavour
0715   if (flavourParton.isNull()) {
0716     if (hardestParton.isNonnull())
0717       partonFlavour = hardestParton->pdgId();
0718   } else
0719     partonFlavour = flavourParton->pdgId();
0720
0721   // if enabled, check for conflicts between hadron- and parton-based flavours and give priority to the hadron-based flavour
0722   if (hadronFlavourHasPriority_) {
0723     if (hadronFlavour == 0 && (std::abs(partonFlavour) == 4 || std::abs(partonFlavour) == 5))
0724       partonFlavour = (hardestLightParton.isNonnull() ? hardestLightParton->pdgId() : 0);
0725     else if (hadronFlavour != 0 && std::abs(partonFlavour) != hadronFlavour)
0726       partonFlavour = hadronFlavour;
0727   }
0728 }
0729
0730 // ------------ method that assigns clustered particles to subjets ------------
0731 void JetFlavourClustering::assignToSubjets(const reco::GenParticleRefVector& clusteredParticles,
0732                                            const edm::Handle<edm::View<reco::Jet>>& subjets,
0733                                            const std::vector<int>& subjetIndices,
0734                                            std::vector<reco::GenParticleRefVector>& assignedParticles) {
0735   // loop over clustered particles and assign them to different subjets based on smallest dR
0736   for (reco::GenParticleRefVector::const_iterator it = clusteredParticles.begin(); it != clusteredParticles.end();
0737        ++it) {
0738     std::vector<double> dR2toSubjets;
0739
0740     for (size_t sj = 0; sj < subjetIndices.size(); ++sj)
0741       dR2toSubjets.push_back(reco::deltaR2((*it)->rapidity(),
0742                                            (*it)->phi(),
0743                                            subjets->at(subjetIndices.at(sj)).rapidity(),
0744                                            subjets->at(subjetIndices.at(sj)).phi()));
0745
0746     // find the closest subjet
0747     int closestSubjetIdx =
0748         std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
0749
0750     assignedParticles.at(closestSubjetIdx).push_back(*it);
0751   }
0752 }
0753
0754 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0755 void JetFlavourClustering::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0756   //The following says we do not know what parameters are allowed so do no validation
0757   // Please change this to state exactly what you do use, even if it is no parameters
0758   edm::ParameterSetDescription desc;
0759   desc.setUnknown();
0760   descriptions.addDefault(desc);
0761 }
0762
0763 //define this as a plug-in
0764 DEFINE_FWK_MODULE(JetFlavourClustering);