JetMCAlgos/plugins/GenHFHadronMatcher.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    GenHFHadronMatcher
0004 // Class:      GenHFHadronMatcher
0005 //
0006
0007 /**\class GenHFHadronMatcher GenHFHadronMatcher.cc
0008 * @brief Finds the origin of each heavy flavour hadron and associated jets to it
0009 *
0010 * Starting from each consituent of each jet, tracks back in chain to find heavy flavour hadrons.
0011 * From each hadron traces back until finds the b quark and its mother.
0012 * For each hadron identifies the jet to which it was injected as a ghost hadron.
0013 *
0014 * The description of the run-time parameters can be found at fillDescriptions()
0015 *
0016 * The description of the products can be found at GenHFHadronMatcher()
0017 */
0018
0019 // Original Author:  Nazar Bartosik,DESY
0020
0021 // system include files
0022 #include <memory>
0023 #include <utility>
0024 #include <vector>
0025 #include <algorithm>
0026
0027 // user include files
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/global/EDProducer.h"
0030
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037 #include "FWCore/Framework/interface/ESHandle.h"
0038 #include "FWCore/Framework/interface/EventSetup.h"
0039 #include "FWCore/Utilities/interface/EDPutToken.h"
0040
0041 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0042
0043 #include "DataFormats/JetReco/interface/GenJet.h"
0044 #include "DataFormats/JetReco/interface/Jet.h"
0045 #include "DataFormats/Math/interface/deltaR.h"
0046
0047 #include "DataFormats/JetMatching/interface/JetFlavourInfoMatching.h"
0048
0049 //
0050 // class declaration
0051 //
0052
0053 class GenHFHadronMatcher : public edm::global::EDProducer<> {
0054 public:
0055   explicit GenHFHadronMatcher(const edm::ParameterSet &);
0056
0057   ~GenHFHadronMatcher() override = default;
0058
0059   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0060
0061 private:
0062   void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0063
0064   std::vector<int> findHadronJets(const reco::GenParticleCollection *genParticles,
0065                                   const reco::JetFlavourInfoMatchingCollection *jetFlavourInfos,
0066                                   std::vector<int> &hadIndex,
0067                                   std::vector<reco::GenParticle> &hadMothersGenPart,
0068                                   std::vector<std::vector<int>> &hadMothersIndices,
0069                                   std::vector<int> &hadLeptonIndex,
0070                                   std::vector<int> &hadLeptonHadIndex,
0071                                   std::vector<int> &hadLeptonViaTau,
0072                                   std::vector<int> &hadFlavour,
0073                                   std::vector<int> &hadFromTopWeakDecay,
0074                                   std::vector<int> &hadBHadronId) const;
0075
0076   int analyzeMothers(const reco::Candidate *thisParticle,
0077                      int &topDaughterQId,
0078                      int &topBarDaughterQId,
0079                      std::vector<const reco::Candidate *> &hadMothers,
0080                      std::vector<std::vector<int>> &hadMothersIndices,
0081                      std::set<const reco::Candidate *> &analyzedParticles,
0082                      const int prevPartIndex) const;
0083   bool putMotherIndex(std::vector<std::vector<int>> &hadMothersIndices, int partIndex, int mothIndex) const;
0084   bool isHadron(const int flavour, const reco::Candidate *thisParticle) const;
0085   bool isHadronPdgId(const int flavour, const int pdgId) const;
0086   bool isMesonPdgId(const int flavour, const int pdgId) const;
0087   bool isBaryonPdgId(const int flavour, const int pdgId) const;
0088   int flavourSign(const int pdgId) const;
0089   bool hasHadronDaughter(const int flavour, const reco::Candidate *thisParticle) const;
0090   int idInList(std::vector<const reco::Candidate *> particleList, const reco::Candidate *particle) const;
0091   int idInList(std::vector<int> list, const int value) const;
0092   int findInMothers(int idx,
0093                     std::vector<int> &mothChains,
0094                     const std::vector<std::vector<int>> &hadMothersIndices,
0095                     const std::vector<reco::GenParticle> &hadMothers,
0096                     int status,
0097                     int pdgId,
0098                     bool pdgAbs,
0099                     int stopId,
0100                     int firstLast,
0101                     bool verbose) const;
0102   bool isNeutralPdg(int pdgId) const;
0103
0104   bool checkForLoop(std::vector<const reco::Candidate *> &particleChain, const reco::Candidate *particle) const;
0105   std::string getParticleName(int id) const;
0106
0107   bool fixExtraSameFlavours(const unsigned int hadId,
0108                             const std::vector<int> &hadIndices,
0109                             const std::vector<reco::GenParticle> &hadMothers,
0110                             const std::vector<std::vector<int>> &hadMothersIndices,
0111                             const std::vector<int> &isFromTopWeakDecay,
0112                             const std::vector<std::vector<int>> &LastQuarkIds,
0113                             const std::vector<std::vector<int>> &LastQuarkMotherIds,
0114                             std::vector<int> &lastQuarkIndices,
0115                             std::vector<int> &hadronFlavour,
0116                             std::set<int> &checkedHadronIds,
0117                             const int lastQuarkIndex) const;
0118
0119   // ----------member data ---------------------------
0120   const edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0121   const edm::EDGetTokenT<reco::JetFlavourInfoMatchingCollection> jetFlavourInfosToken_;
0122   const int flavour_;
0123   const bool noBBbarResonances_;
0124   const bool onlyJetClusteredHadrons_;
0125
0126   const std::string flavourStr_;
0127   const edm::EDPutTokenT<std::vector<reco::GenParticle>> plusMothersToken_;
0128   const edm::EDPutTokenT<std::vector<std::vector<int>>> plusMothersIndicesToken_;
0129   const edm::EDPutTokenT<std::vector<int>> indexToken_;
0130   const edm::EDPutTokenT<std::vector<int>> flavourToken_;
0131   const edm::EDPutTokenT<std::vector<int>> jetIndexToken_;
0132   const edm::EDPutTokenT<std::vector<int>> leptonIndexToken_;
0133   const edm::EDPutTokenT<std::vector<int>> leptonHadronIndexToken_;
0134   const edm::EDPutTokenT<std::vector<int>> leptonViaTauToken_;
0135   const edm::EDPutTokenT<std::vector<int>> fromTopWeakDecayToken_;
0136   const edm::EDPutTokenT<std::vector<int>> bHadronIdToken_;
0137 };
0138
0139 //
0140 // constants, enums and typedefs
0141 //
0142
0143 //
0144 // static data member definitions
0145 //
0146
0147 //
0148 // constructors and destructor
0149 //
0150 /**
0151 * @brief constructor initialising producer products and config parameters
0152 *
0153 * @warning Definition of b-hadron and anti-b-hadron: The term b-hadron and anti-b-hadron is in reference to the quark content and <b>not</b> to distinguish particles from anti-particles.
0154 * Here a b-hadron contains a b-quark and an anti-b-hadron contains an anti-b-quark.
0155 * For mesons this means an inversion with respect to the PDG definition, as mesons actually contain anti-b-quarks and anti-mesons contain b-quarks.
0156 *
0157 */
0158 namespace {
0159   std::string flavourName(int flavour) {
0160     if (flavour == 5) {
0161       return "B";
0162     } else if (flavour == 4) {
0163       return "C";
0164     }
0165     edm::LogError("GenHFHadronMatcher") << "Flavour option must be 4 (c-jet) or 5 (b-jet), but is: " << flavour
0166                                         << ". Correct this!";
0167     return std::string();
0168   }
0169 }  // namespace
0170
0171 GenHFHadronMatcher::GenHFHadronMatcher(const edm::ParameterSet &cfg)
0172     : genParticlesToken_(consumes<reco::GenParticleCollection>(cfg.getParameter<edm::InputTag>("genParticles"))),
0173       jetFlavourInfosToken_(
0174           consumes<reco::JetFlavourInfoMatchingCollection>(cfg.getParameter<edm::InputTag>("jetFlavourInfos"))),
0175       flavour_{std::abs(cfg.getParameter<int>("flavour"))},
0176       noBBbarResonances_{cfg.getParameter<bool>("noBBbarResonances")},
0177       onlyJetClusteredHadrons_{cfg.getParameter<bool>("onlyJetClusteredHadrons")},
0178       flavourStr_{flavourName(flavour_)},
0179       plusMothersToken_{produces<std::vector<reco::GenParticle>>(
0180           "gen" + flavourStr_ +
0181           "HadPlusMothers")},  // All mothers in all decay chains above any hadron of specified flavour
0182       plusMothersIndicesToken_{produces<std::vector<std::vector<int>>>(
0183           "gen" + flavourStr_ + "HadPlusMothersIndices")},  // Indices of mothers of each hadMother
0184       indexToken_{
0185           produces<std::vector<int>>("gen" + flavourStr_ + "HadIndex")},  // Index of hadron in the vector of hadMothers
0186       flavourToken_{produces<std::vector<int>>(
0187           "gen" + flavourStr_ +
0188           "HadFlavour")},  // PdgId of the first non-b(c) quark mother with sign corresponding to hadron charge
0189       jetIndexToken_{produces<std::vector<int>>(
0190           "gen" + flavourStr_ + "HadJetIndex")},  // Index of genJet matched to each hadron by jet clustering algorithm
0191       leptonIndexToken_{produces<std::vector<int>>(
0192           "gen" + flavourStr_ +
0193           "HadLeptonIndex")},  // Index of lepton found among the hadron decay products in the list of mothers
0194       leptonHadronIndexToken_{produces<std::vector<int>>(
0195           "gen" + flavourStr_ + "HadLeptonHadronIndex")},  // Index of hadron the lepton is associated to
0196       leptonViaTauToken_{produces<std::vector<int>>(
0197           "gen" + flavourStr_ + "HadLeptonViaTau")},  // Whether lepton comes directly from hadron or via tau decay
0198       fromTopWeakDecayToken_{produces<std::vector<int>>(
0199           "gen" + flavourStr_ +
0200           "HadFromTopWeakDecay")},  // Tells whether the hadron appears in the chain after top decay
0201       bHadronIdToken_{produces<std::vector<int>>("gen" + flavourStr_ + "HadBHadronId")}
0202 // Index of a b-hadron which the current hadron comes from (for c-hadrons)
0203 {
0204   // Hadron matching products
0205 }
0206
0207 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0208 /**
0209 * @brief description of the run-time parameters
0210 *
0211 */
0212 void GenHFHadronMatcher::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0213   edm::ParameterSetDescription desc;
0214   desc.add<edm::InputTag>("genParticles")
0215       ->setComment("Collection of GenParticle objects which contains all particles produced in the event");
0216   desc.add<edm::InputTag>("jetFlavourInfos")
0217       ->setComment(
0218           "Output from the JetFlavour tool. Contains information about partons/hadrons/leptons associated to jets");
0219   desc.add<bool>("noBBbarResonances", true)->setComment("Whether resonances should not be treated as hadrons");
0220   desc.add<bool>("onlyJetClusteredHadrons", false)
0221       ->setComment("Whether only hadrons that are matched to jets should be analysed. Runs x1000 faster in Sherpa");
0222   desc.add<int>("flavour", 5)->setComment("Flavour of weakly decaying hadron that should be analysed (4-c, 5-b)");
0223   descriptions.add("matchGenHFHadron", desc);
0224 }
0225
0226 //
0227 // member functions
0228 //
0229
0230 // ------------ method called to produce the data  ------------
0231 void GenHFHadronMatcher::produce(edm::StreamID, edm::Event &evt, const edm::EventSetup &setup) const {
0232   using namespace edm;
0233
0234   edm::Handle<reco::GenParticleCollection> genParticles;
0235   evt.getByToken(genParticlesToken_, genParticles);
0236
0237   edm::Handle<reco::JetFlavourInfoMatchingCollection> jetFlavourInfos;
0238   evt.getByToken(jetFlavourInfosToken_, jetFlavourInfos);
0239
0240   // Defining hadron matching variables
0241   std::vector<reco::GenParticle> hadMothers;
0242   std::vector<std::vector<int>> hadMothersIndices;
0243   std::vector<int> hadIndex;
0244   std::vector<int> hadFlavour;
0245   std::vector<int> hadJetIndex;
0246   std::vector<int> hadLeptonIndex;
0247   std::vector<int> hadLeptonHadIndex;
0248   std::vector<int> hadLeptonViaTau;
0249   std::vector<int> hadFromTopWeakDecay;
0250   std::vector<int> hadBHadronId;
0251
0252   hadJetIndex = findHadronJets(genParticles.product(),
0253                                jetFlavourInfos.product(),
0254                                hadIndex,
0255                                hadMothers,
0256                                hadMothersIndices,
0257                                hadLeptonIndex,
0258                                hadLeptonHadIndex,
0259                                hadLeptonViaTau,
0260                                hadFlavour,
0261                                hadFromTopWeakDecay,
0262                                hadBHadronId);
0263
0264   // Putting products to the event
0265   evt.emplace(plusMothersToken_, std::move(hadMothers));
0266   evt.emplace(plusMothersIndicesToken_, std::move(hadMothersIndices));
0267   evt.emplace(indexToken_, std::move(hadIndex));
0268   evt.emplace(flavourToken_, std::move(hadFlavour));
0269   evt.emplace(jetIndexToken_, std::move(hadJetIndex));
0270   evt.emplace(leptonIndexToken_, std::move(hadLeptonIndex));
0271   evt.emplace(leptonHadronIndexToken_, std::move(hadLeptonHadIndex));
0272   evt.emplace(leptonViaTauToken_, std::move(hadLeptonViaTau));
0273   evt.emplace(fromTopWeakDecayToken_, std::move(hadFromTopWeakDecay));
0274   evt.emplace(bHadronIdToken_, std::move(hadBHadronId));
0275 }
0276
0277 /**
0278 * @brief identify the jets that contain b-hadrons
0279 *
0280 * For each hadron all mothers from all levels and chains are analysed to find the quark or gluon from which the hadron has originated.
0281 * This is done by searching through the generator particle decay tree starting from the hadron, performing checks for flavour and kinematic consistency.
0282 * The hadrons that are not last in the decay chain (don't decay weakly) are skipped.
0283 * Additionally for each hadron it is checked whether it comes from the top weak decay branch and whether it comes from some other b-hadron decay
0284 *
0285 * b-bbar (c-cbar) resonances can either be considered as hadrons or not, depending on the configuration.
0286 *
0287 * @param[out] hadIndex vector of indices of found hadrons in hadMothers
0288 * @param[out] hadMothers vector of all mothers at all levels of each found hadron
0289 * @param[out] hadMothersIndices connection between each particle from hadMothers and its mothers
0290 * @param[out] hadLeptonIndex index of lepton among the hadMothers
0291 * @param[out] hadLeptonHadIndex index of hadron associated to each lepton
0292 * @param[out] hadLeptonViaTau whether lepton is from direct hadron->lepton or hadron->tau->lepton
0293 * @param[out] hadFlavour flavour of each found hadron
0294 * @param[out] hadFromTopWeakDecay whether each hadron contains the top quark in its decay chain [works only for B-Hadrons]
0295 * @param[out] hadBHadronId for each hadron - index of the ancestor b-hadron [-1 if hadron doesn't come from another b-hdaron]
0296 *
0297 * @returns vector of jet indices that were matched to each hadron [by the jet clustering algorithm]
0298 */
0299 std::vector<int> GenHFHadronMatcher::findHadronJets(const reco::GenParticleCollection *genParticles,
0300                                                     const reco::JetFlavourInfoMatchingCollection *jetFlavourInfos,
0301                                                     std::vector<int> &hadIndex,
0302                                                     std::vector<reco::GenParticle> &hadMothers,
0303                                                     std::vector<std::vector<int>> &hadMothersIndices,
0304                                                     std::vector<int> &hadLeptonIndex,
0305                                                     std::vector<int> &hadLeptonHadIndex,
0306                                                     std::vector<int> &hadLeptonViaTau,
0307                                                     std::vector<int> &hadFlavour,
0308                                                     std::vector<int> &hadFromTopWeakDecay,
0309                                                     std::vector<int> &hadBHadronId) const {
0310   std::vector<int> hadJetIndex;
0311   std::vector<const reco::Candidate *> hadMothersCand;
0312
0313   int topDaughterQId = -1;
0314   int topBarDaughterQId = -1;
0315
0316   // Looping over all jets to get hadrons associated to them
0317   for (reco::JetFlavourInfoMatchingCollection::const_iterator i_info = jetFlavourInfos->begin();
0318        i_info != jetFlavourInfos->end();
0319        ++i_info) {
0320     reco::JetFlavourInfo jetInfo = i_info->second;
0321     const int jetIndex = i_info - jetFlavourInfos->begin();
0322     // Looping over each hadron associated with the jet and finding its origin
0323     const reco::GenParticleRefVector &hadronsInJet = flavour_ == 5   ? jetInfo.getbHadrons()
0324                                                      : flavour_ == 4 ? jetInfo.getcHadrons()
0325                                                                      : reco::GenParticleRefVector();
0326     for (reco::GenParticleRefVector::const_iterator hadron = hadronsInJet.begin(); hadron != hadronsInJet.end();
0327          ++hadron) {
0328       // Check that the hadron satisfies criteria configured in the module
0329       if (!isHadron(flavour_, (&**hadron)))
0330         continue;
0331       if (hasHadronDaughter(flavour_, (reco::Candidate *)(&**hadron)))
0332         continue;
0333       // Scanning the chain starting from the hadron
0334       std::set<const reco::Candidate *> analyzedParticles;
0335       int hadronIndex = analyzeMothers((reco::Candidate *)(&**hadron),
0336                                        topDaughterQId,
0337                                        topBarDaughterQId,
0338                                        hadMothersCand,
0339                                        hadMothersIndices,
0340                                        analyzedParticles,
0341                                        -1);
0342       // Storing the index of the hadron to the list
0343       hadIndex.push_back(hadronIndex);
0344       hadJetIndex.push_back(jetIndex);  // Putting jet index to the result list
0345     }
0346   }  // End of loop over jets
0347
0348   // Access all hadrons which are not associated with jets, if requested
0349   if (!onlyJetClusteredHadrons_) {
0350     for (reco::GenParticleCollection::const_iterator i_particle = genParticles->begin();
0351          i_particle != genParticles->end();
0352          ++i_particle) {
0353       const reco::GenParticle *thisParticle = &*i_particle;
0354       if (!isHadron(flavour_, thisParticle))
0355         continue;
0356       // Skipping the hadron if it was already found directly from jets
0357       if (std::find(hadMothersCand.begin(), hadMothersCand.end(), thisParticle) != hadMothersCand.end())
0358         continue;
0359
0360       // Scanning the chain starting from the hadron
0361       std::set<const reco::Candidate *> analyzedParticles;
0362       int hadronIndex = analyzeMothers(
0363           thisParticle, topDaughterQId, topBarDaughterQId, hadMothersCand, hadMothersIndices, analyzedParticles, -1);
0364       // Storing the index of the hadron to the list
0365       hadIndex.push_back(hadronIndex);
0366       hadJetIndex.push_back(-1);  // Jet index undefined
0367     }
0368   }
0369
0370   // Transfering Candidates to the list of processed particles for further analysis
0371   for (int i = 0; i < (int)hadMothersCand.size(); i++) {
0372     const reco::GenParticle *particle = dynamic_cast<const reco::GenParticle *>(hadMothersCand.at(i));
0373     hadMothers.push_back(*particle);
0374   }
0375
0376   // Adding leptons from hadron decays
0377   for (reco::GenParticleCollection::const_iterator i_particle = genParticles->begin();
0378        i_particle != genParticles->end();
0379        ++i_particle) {
0380     const reco::GenParticle &lepton = *i_particle;
0381     const int pdg_abs = lepton.pdgId();
0382     // Skipping if not a lepton: e/mu
0383     if (pdg_abs != 11 && pdg_abs != 13)
0384       continue;
0385     bool leptonViaTau = false;
0386     const reco::Candidate *leptonMother = lepton.mother();
0387     if (!leptonMother)
0388       continue;
0389     // Taking next mother if direct mother is a tau
0390     if (std::abs(leptonMother->pdgId()) == 15) {
0391       leptonViaTau = true;
0392       leptonMother = leptonMother->mother();
0393     }
0394     // Skipping this lepton if its mother is not a proper hadron
0395     if (leptonMother == nullptr or !isHadron(flavour_, leptonMother))
0396       continue;
0397     // Finding the index of this hadron in the list of analysed particles
0398     size_t leptonHadronParticleIndex =
0399         std::find(hadMothersCand.begin(), hadMothersCand.end(), leptonMother) - hadMothersCand.begin();
0400     if (leptonHadronParticleIndex >= hadMothersCand.size())
0401       continue;
0402     // Finding the actual hadron index among those that were found
0403     size_t leptonHadronIndex =
0404         std::find(hadIndex.begin(), hadIndex.end(), leptonHadronParticleIndex) - hadIndex.begin();
0405     if (leptonHadronIndex >= hadIndex.size())
0406       continue;
0407     // Putting the lepton, its index and hadron index to the corresponding lists
0408     hadMothers.push_back(lepton);
0409     const int leptonIndex = hadMothersCand.size() - 1;
0410     hadLeptonIndex.push_back(leptonIndex);
0411     hadLeptonViaTau.push_back(leptonViaTau);
0412     hadLeptonHadIndex.push_back(leptonHadronIndex);
0413   }
0414
0415   // Checking mothers of hadrons in order to assign flags (where the hadron comes from)
0416   unsigned int nHad = hadIndex.size();
0417
0418   std::vector<std::vector<int>> LastQuarkMotherIds;
0419   std::vector<std::vector<int>> LastQuarkIds;
0420   std::vector<int> lastQuarkIndices(nHad, -1);
0421
0422   // Looping over all hadrons
0423   for (unsigned int hadNum = 0; hadNum < nHad; hadNum++) {
0424     unsigned int hadIdx = hadIndex.at(hadNum);  // Index of hadron in the hadMothers
0425
0426     std::vector<int> FirstQuarkId;
0427     std::vector<int> LastQuarkId;
0428     std::vector<int> LastQuarkMotherId;
0429
0430     const int hadronFlavourSign = flavourSign(hadMothers.at(hadIdx).pdgId());
0431
0432     // Searching only first quark in the chain with the same flavour as hadron
0433     if (hadronFlavourSign != 0) {
0434       findInMothers(
0435           hadIdx, FirstQuarkId, hadMothersIndices, hadMothers, 0, hadronFlavourSign * flavour_, false, -1, 1, false);
0436     }
0437     // Searching for quarks with both flavours since it is a bb/cc resonance
0438     else {
0439       findInMothers(hadIdx, FirstQuarkId, hadMothersIndices, hadMothers, 0, flavour_, false, -1, 1, false);
0440       findInMothers(hadIdx, FirstQuarkId, hadMothersIndices, hadMothers, 0, -1 * flavour_, false, -1, 1, false);
0441     }
0442
0443     // Finding last quark for each first quark
0444     for (unsigned int qId = 0; qId < FirstQuarkId.size(); qId++) {
0445       // Identifying the flavour of the first quark to find the last quark of the same flavour
0446       const int quarkFlavourSign = flavourSign(hadMothers.at(FirstQuarkId.at(qId)).pdgId());
0447       // Finding last quark of the hadron starting from the first quark
0448       findInMothers(FirstQuarkId.at(qId),
0449                     LastQuarkId,
0450                     hadMothersIndices,
0451                     hadMothers,
0452                     0,
0453                     quarkFlavourSign * flavour_,
0454                     false,
0455                     -1,
0456                     2,
0457                     false);
0458     }  // End of loop over all first quarks of the hadron
0459
0460     // Setting initial flavour of the hadron
0461     int hadronFlavour = 0;
0462
0463     // Initialising pairs of last quark index and its distance from the hadron (to sort by distance)
0464     std::vector<std::pair<double, int>> lastQuark_dR_id_pairs;
0465
0466     // Finding the closest quark in dR
0467     for (unsigned int qId = 0; qId < LastQuarkId.size(); qId++) {
0468       int qIdx = LastQuarkId.at(qId);
0469       // Calculating the dR between hadron and quark
0470       float dR = deltaR(hadMothers.at(hadIdx).eta(),
0471                         hadMothers.at(hadIdx).phi(),
0472                         hadMothers.at(qIdx).eta(),
0473                         hadMothers.at(qIdx).phi());
0474
0475       std::pair<double, int> dR_hadId_pair(dR, qIdx);
0476       lastQuark_dR_id_pairs.push_back(dR_hadId_pair);
0477     }  // End of loop over all last quarks of the hadron
0478
0479     std::sort(lastQuark_dR_id_pairs.begin(), lastQuark_dR_id_pairs.end());
0480
0481     if (lastQuark_dR_id_pairs.size() > 1) {
0482       double dRratio =
0483           (lastQuark_dR_id_pairs.at(1).first - lastQuark_dR_id_pairs.at(0).first) / lastQuark_dR_id_pairs.at(1).first;
0484       int qIdx_closest = lastQuark_dR_id_pairs.at(0).second;
0485       LastQuarkId.clear();
0486       if (dRratio > 0.5)
0487         LastQuarkId.push_back(qIdx_closest);
0488       else
0489         for (std::pair<double, int> qIdDrPair : lastQuark_dR_id_pairs)
0490           LastQuarkId.push_back(qIdDrPair.second);
0491     }
0492     for (int qIdx : LastQuarkId) {
0493       int qmIdx = hadMothersIndices.at(qIdx).at(0);
0494       LastQuarkMotherId.push_back(qmIdx);
0495     }
0496
0497     if ((int)LastQuarkId.size() > 0)
0498       lastQuarkIndices.at(hadNum) = 0;  // Setting the first quark in array as a candidate if it exists
0499
0500     LastQuarkIds.push_back(LastQuarkId);
0501
0502     LastQuarkMotherIds.push_back(LastQuarkMotherId);
0503
0504     if (LastQuarkMotherId.empty()) {
0505       hadronFlavour = 0;
0506     } else {
0507       int qIdx = LastQuarkId.at(lastQuarkIndices.at(hadNum));
0508       int qFlav = flavourSign(hadMothers.at(qIdx).pdgId());
0509       hadronFlavour = qFlav * std::abs(hadMothers.at(LastQuarkMotherId.at(lastQuarkIndices.at(hadNum))).pdgId());
0510     }
0511     hadFlavour.push_back(hadronFlavour);  // Adding hadron flavour to the list of flavours
0512
0513     // Checking whether hadron comes from the Top weak decay
0514     int isFromTopWeakDecay = 1;
0515     std::vector<int> checkedParticles;
0516     if (hadFlavour.at(hadNum) != 0) {
0517       int lastQIndex = LastQuarkId.at(lastQuarkIndices.at(hadNum));
0518       bool fromTB = topDaughterQId >= 0 ? findInMothers(lastQIndex,
0519                                                         checkedParticles,
0520                                                         hadMothersIndices,
0521                                                         hadMothers,
0522                                                         -1,
0523                                                         0,
0524                                                         false,
0525                                                         topDaughterQId,
0526                                                         2,
0527                                                         false) >= 0
0528                                         : false;
0529       checkedParticles.clear();
0530       bool fromTbarB = topBarDaughterQId >= 0 ? findInMothers(lastQIndex,
0531                                                               checkedParticles,
0532                                                               hadMothersIndices,
0533                                                               hadMothers,
0534                                                               -1,
0535                                                               0,
0536                                                               false,
0537                                                               topBarDaughterQId,
0538                                                               2,
0539                                                               false) >= 0
0540                                               : false;
0541       checkedParticles.clear();
0542       if (!fromTB && !fromTbarB) {
0543         isFromTopWeakDecay = 0;
0544       }
0545     } else
0546       isFromTopWeakDecay = 2;
0547     hadFromTopWeakDecay.push_back(isFromTopWeakDecay);
0548     int bHadronMotherId =
0549         findInMothers(hadIdx, checkedParticles, hadMothersIndices, hadMothers, 0, 555555, true, -1, 1, false);
0550     hadBHadronId.push_back(bHadronMotherId);
0551
0552     if (!LastQuarkMotherId.empty()) {
0553       std::set<int> checkedHadronIds;
0554       fixExtraSameFlavours(hadNum,
0555                            hadIndex,
0556                            hadMothers,
0557                            hadMothersIndices,
0558                            hadFromTopWeakDecay,
0559                            LastQuarkIds,
0560                            LastQuarkMotherIds,
0561                            lastQuarkIndices,
0562                            hadFlavour,
0563                            checkedHadronIds,
0564                            0);
0565     }
0566
0567   }  // End of loop over all hadrons
0568
0569   return hadJetIndex;
0570 }
0571
0572 /**
0573 * @brief Check if the cpecified particle is already in the list of particles
0574 *
0575 * @param[in] particleList list of particles to be checked
0576 * @param[in] particle particle that should be checked
0577 *
0578 * @returns the index of the particle in the list [-1 if particle not found]
0579 */
0580 int GenHFHadronMatcher::idInList(std::vector<const reco::Candidate *> particleList,
0581                                  const reco::Candidate *particle) const {
0582   const unsigned int position = std::find(particleList.begin(), particleList.end(), particle) - particleList.begin();
0583   if (position >= particleList.size())
0584     return -1;
0585
0586   return position;
0587 }
0588
0589 int GenHFHadronMatcher::idInList(std::vector<int> list, const int value) const {
0590   const unsigned int position = std::find(list.begin(), list.end(), value) - list.begin();
0591   if (position >= list.size())
0592     return -1;
0593
0594   return position;
0595 }
0596
0597 /**
0598 * @brief Check the pdgId of a given particle if it is a hadron
0599 *
0600 * @param[in] flavour flavour of a hadron that is being searched (5-B, 4-C)
0601 * @param[in] thisParticle a particle that is to be analysed
0602 *
0603 * @returns whether the particle is a hadron of specified flavour
0604 */
0605 bool GenHFHadronMatcher::isHadron(const int flavour, const reco::Candidate *thisParticle) const {
0606   return isHadronPdgId(flavour, thisParticle->pdgId());
0607 }
0608
0609 /**
0610 * @brief Check the pdgId if it represents a hadron of particular flavour
0611 *
0612 * @param[in] flavour flavour of a hadron that is being searched (5-B, 4-C)
0613 * @param[in] pdgId pdgId to be checked
0614 *
0615 * @returns true if the pdgId represents a hadron of specified flavour
0616 */
0617 bool GenHFHadronMatcher::isHadronPdgId(const int flavour, const int pdgId) const {
0618   if (isBaryonPdgId(flavour, pdgId) || isMesonPdgId(flavour, pdgId))
0619     return true;
0620
0621   return false;
0622 }
0623
0624 /**
0625 * @brief Check the pdgId if it represents a meson of particular flavour
0626 *
0627 * @param[in] flavour flavour of a hadron that is being searched (5-B, 4-C)
0628 * @param[in] pdgId pdgId to be checked
0629 *
0630 * @returns true if the pdgId represents a meson of specified flavour
0631 */
0632 bool GenHFHadronMatcher::isMesonPdgId(const int flavour, const int pdgId) const {
0633   const int flavour_abs = std::abs(flavour);
0634   if (flavour_abs != 5 && flavour_abs != 4)
0635     return false;
0636   const int pdgId_abs = std::abs(pdgId);
0637
0638   if (pdgId_abs / 100 % 10 != flavour_abs)
0639     return false;
0640   // Excluding baryons
0641   if (pdgId_abs / 1000 == flavour_abs)
0642     return false;
0643   // Excluding bb/cc resonances if required
0644   if (noBBbarResonances_ && pdgId_abs / 10 % 100 == 11 * flavour_abs)
0645     return false;
0646
0647   return true;
0648 }
0649
0650 /**
0651 * @brief Check the pdgId if it represents a baryon of particular flavour
0652 *
0653 * @param[in] flavour flavour of a hadron that is being searched (5-B, 4-C)
0654 * @param[in] pdgId pdgId to be checked
0655 *
0656 * @returns true if the pdgId represents a baryon of specified flavour
0657 */
0658 bool GenHFHadronMatcher::isBaryonPdgId(const int flavour, const int pdgId) const {
0659   const int flavour_abs = std::abs(flavour);
0660   if (flavour_abs != 5 && flavour_abs != 4)
0661     return false;
0662   const int pdgId_abs = std::abs(pdgId);
0663
0664   if (pdgId_abs / 1000 != flavour_abs)
0665     return false;
0666
0667   return true;
0668 }
0669
0670 /**
0671 * @brief Sign of the flavour (matter/antimatter)
0672 *
0673 * @param[in] pdgId pdgId to be checked
0674 *
0675 * @returns +1/-1/0  matter/antimatter/undefined
0676 */
0677 int GenHFHadronMatcher::flavourSign(const int pdgId) const {
0678   int flavourSign = pdgId / std::abs(pdgId);
0679   // B mesons have opposite sign
0680   if (isMesonPdgId(5, pdgId))
0681     flavourSign *= -1;
0682   // Returning 0 for bb/cc resonances
0683   if (pdgId % 1000 / 10 / 11 > 0)
0684     flavourSign = 0;
0685
0686   return flavourSign;
0687 }
0688
0689 /**
0690 * @brief Check if the particle has bHadron among daughters
0691 *
0692 * @param[in] flavour flavour of a hadron that is being searched (5-B, 4-C)
0693 * @param[in] thisParticle a particle that is to be analysed
0694 *
0695 * @returns whether the particle has a hadron among its daughters
0696 */
0697 bool GenHFHadronMatcher::hasHadronDaughter(const int flavour, const reco::Candidate *thisParticle) const {
0698   // Looping through daughters of the particle
0699   bool hasDaughter = false;
0700   for (int k = 0; k < (int)thisParticle->numberOfDaughters(); k++) {
0701     if (!isHadron(flavour, thisParticle->daughter(k))) {
0702       continue;
0703     }
0704     hasDaughter = true;
0705     break;
0706   }
0707
0708   return hasDaughter;
0709 }
0710
0711 /**
0712 * @brief do a recursive search for the mother particles until the b-quark is found or the absolute mother is found
0713 *
0714 * the treatment of b-bar resonances depends on the global parameter noBBbarResonances_
0715 *
0716 * @param[in] thisParticle current particle from which starts the search of the hadron and all its mothers up to proton
0717 * @param[out] hadron the last hadron in the decay chain [that decays weekly]
0718 * @param[out] lepton lepton found in the current decay chain
0719 * @param[out] topDaughterQId Id of the top quark daughter b(c) quark
0720 * @param[out] topBarDaughterQId Id of the antitop quark daughter b(c) quark
0721 * @param[out] hadMothers list of all processed particles ending with proton
0722 * @param[out] hadMothersIndices list of i-vectors containing j-indices representing particles that are mothers of each i-particle from hadMothers
0723 * @param[out] analyzedParticles list of particles analysed in the chain [used for infinite loop detection]
0724 * @param[out] prevPartIndex index of the previous particle in the current chain [used for infinite loop detection]
0725 *
0726 * @returns index of hadron in the hadMothers list [-1 if no hadron found]
0727 */
0728 int GenHFHadronMatcher::analyzeMothers(const reco::Candidate *thisParticle,
0729                                        int &topDaughterQId,
0730                                        int &topBarDaughterQId,
0731                                        std::vector<const reco::Candidate *> &hadMothers,
0732                                        std::vector<std::vector<int>> &hadMothersIndices,
0733                                        std::set<const reco::Candidate *> &analyzedParticles,
0734                                        const int prevPartIndex) const {
0735   // Getting the index of the particle which is a hadron in the first call
0736   int hadronIndex = -1;  // Index of the hadron that is returned by this function
0737   int index = idInList(hadMothers, thisParticle);
0738   if (index < 0) {  // If hadron is not in the list of mothers yet
0739     hadMothers.push_back(thisParticle);
0740     hadronIndex = hadMothers.size() - 1;
0741   } else {  // If hadron is in the list of mothers already
0742     hadronIndex = index;
0743   }
0744
0745   int partIndex = -1;  // Index of particle being checked in the list of mothers
0746   partIndex = idInList(hadMothers, thisParticle);
0747
0748   // Checking whether this particle is already in the chain of analyzed particles in order to identify a loop
0749   bool isLoop = false;
0750   for (unsigned int i = 0; i < analyzedParticles.size(); i++) {
0751     if (analyzedParticles.count(thisParticle) <= 0) {
0752       continue;
0753     }
0754     isLoop = true;
0755     break;
0756   }
0757
0758   // If a loop has been detected
0759   if (isLoop) {
0760     if (prevPartIndex >= 0) {
0761       putMotherIndex(hadMothersIndices, prevPartIndex, -1);  // Setting mother index of previous particle to -1
0762     }
0763     return hadronIndex;  // Stopping further processing of the current chain
0764   }
0765   analyzedParticles.insert(thisParticle);
0766
0767   // Putting the mothers to the list of mothers
0768   for (size_t iMother = 0; iMother < thisParticle->numberOfMothers(); ++iMother) {
0769     const reco::Candidate *mother = thisParticle->mother(iMother);
0770     int mothIndex = idInList(hadMothers, mother);
0771     if (mothIndex == partIndex && partIndex >= 0) {
0772       continue;  // Skipping the mother that is its own daughter
0773     }
0774
0775     // If this mother isn't yet in the list and hadron or lepton is in the list
0776     if (mothIndex < 0) {
0777       hadMothers.push_back(mother);
0778       mothIndex = hadMothers.size() - 1;
0779     }
0780     // If hadron has already been found in current chain and the mother isn't a duplicate of the particle being checked
0781     if (mothIndex != partIndex && partIndex >= 0) {
0782       putMotherIndex(hadMothersIndices, partIndex, mothIndex);  // Putting the index of mother for current particle
0783     }
0784     analyzeMothers(
0785         mother, topDaughterQId, topBarDaughterQId, hadMothers, hadMothersIndices, analyzedParticles, partIndex);
0786     // Setting the id of the particle which is a quark from the top decay
0787     if (std::abs(mother->pdgId()) == 6) {
0788       int &bId = mother->pdgId() < 0 ? topBarDaughterQId : topDaughterQId;
0789       int thisFlav = std::abs(thisParticle->pdgId());
0790       if (bId < 0) {
0791         if (thisFlav <= 5)
0792           bId = partIndex;
0793       } else {
0794         int bIdFlav = std::abs(hadMothers.at(bId)->pdgId());
0795         if (bIdFlav != 5 && thisFlav == 5)
0796           bId = partIndex;
0797         else if (thisFlav == 5 && thisParticle->pt() > hadMothers.at(bId)->pt())
0798           bId = partIndex;
0799       }  // If daughter quark of the top not found yet
0800     }  // If the mother is a top quark and hadron has been found
0801   }  // End of loop over mothers
0802
0803   analyzedParticles.erase(thisParticle);  // Removing current particle from the current chain that is being analyzed
0804
0805   if (partIndex < 0) {
0806     return hadronIndex;  // Safety check
0807   }
0808
0809   // Adding -1 to the list of mother indices for current particle if it has no mothers (for consistency between numbering of indices and mothers)
0810   if ((int)thisParticle->numberOfMothers() <= 0) {
0811     putMotherIndex(hadMothersIndices, partIndex, -1);
0812   }
0813
0814   return hadronIndex;
0815 }
0816
0817 /**
0818 * @brief puts mother index to the list of mothers of particle, if it isn't there already
0819 *
0820 * @param[in] hadMothersIndices vector of indices of mothers for each particle
0821 * @param[in] partIndex index of the particle for which the mother index should be stored
0822 * @param[in] mothIndex index of mother that should be stored for current particle
0823 *
0824 * @returns whether the particle index was alreade in the list
0825 */
0826 bool GenHFHadronMatcher::putMotherIndex(std::vector<std::vector<int>> &hadMothersIndices,
0827                                         int partIndex,
0828                                         int mothIndex) const {
0829   // Putting vector of mothers indices for the given particle
0830   bool inList = false;
0831   if (partIndex < 0) {
0832     return false;
0833   }
0834
0835   while ((int)hadMothersIndices.size() <= partIndex) {  // If there is no list of mothers for current particle yet
0836     std::vector<int> mothersIndices;
0837     hadMothersIndices.push_back(mothersIndices);
0838   }
0839
0840   std::vector<int> *hadMotherIndices = &hadMothersIndices.at(partIndex);
0841   // Removing other mothers if particle must have no mothers
0842   if (mothIndex == -1) {
0843     hadMotherIndices->clear();
0844   } else {
0845     // Checking if current mother is already in the list of theParticle's mothers
0846     for (int k = 0; k < (int)hadMotherIndices->size(); k++) {
0847       if (hadMotherIndices->at(k) != mothIndex && hadMotherIndices->at(k) != -1) {
0848         continue;
0849       }
0850       inList = true;
0851       break;
0852     }
0853   }
0854   // Adding current mother to the list of mothers of this particle
0855   if (!inList) {
0856     hadMotherIndices->push_back(mothIndex);
0857   }
0858
0859   return inList;
0860 }
0861
0862 /**
0863 * @brief helper function to find indices of particles with particular pdgId and status from the list of mothers of a given particle
0864 *
0865 * @param[in] idx index of particle, mothers of which should be searched
0866 * @param[in] mothChains vector of indices where the found mothers are stored
0867 * @param[in] hadMothersIndices list of indices pointing to mothers of each particle from list of mothers
0868 * @param[in] hadMothers vector of all hadron mother particles of all levels
0869 * @param[in] status status of mother that is being looked for
0870 * @param[in] pdgId pdgId of mother that is being looked for [flavour*111111 used to identify hadrons of particular flavour]
0871 * @param[in] pdgAbs whether the sign of pdgId should be taken into account
0872 * @param[in] stopId id of the particle in the hadMothers array after which the checking should stop
0873 * @param[in] firstLast should all(0), first(1) or last(2) occurances of the searched particle be stored
0874 * @param[in] verbose option to print every particle that is processed during the search
0875 *
0876 * @returns index of the found particle in the hadMothers array [-1 if the specified particle not found]
0877 */
0878
0879 int GenHFHadronMatcher::findInMothers(int idx,
0880                                       std::vector<int> &mothChains,
0881                                       const std::vector<std::vector<int>> &hadMothersIndices,
0882                                       const std::vector<reco::GenParticle> &hadMothers,
0883                                       int status,
0884                                       int pdgId,
0885                                       bool pdgAbs = false,
0886                                       int stopId = -1,
0887                                       int firstLast = 0,
0888                                       bool verbose = false) const {
0889   int foundStopId = -1;
0890   int pdg_1 = hadMothers.at(idx).pdgId();
0891
0892   if ((int)hadMothersIndices.size() <= idx) {
0893     if (verbose) {
0894       printf(" Stopping checking particle %d. No mothers are stored.\n", idx);
0895     }
0896     return -1;  // Skipping if no mothers are stored for this particle
0897   }
0898
0899   if (std::abs(hadMothers.at(idx).pdgId()) > 10 && std::abs(hadMothers.at(idx).pdgId()) < 19)
0900     printf("Lepton: %d\n", hadMothers.at(idx).pdgId());
0901
0902   std::vector<int> mothers = hadMothersIndices.at(idx);
0903   unsigned int nMothers = mothers.size();
0904   bool isCorrect = false;  // Whether current particle is what is being searched
0905   if (verbose) {
0906     if (std::abs(hadMothers.at(idx).pdgId()) == 2212) {
0907       printf("Chk:  %d\tpdg: %d\tstatus: %d", idx, hadMothers.at(idx).pdgId(), hadMothers.at(idx).status());
0908     } else {
0909       printf(" Chk:  %d(%d mothers)\tpdg: %d\tstatus: %d\tPt: %.3f\tEta: %.3f",
0910              idx,
0911              nMothers,
0912              hadMothers.at(idx).pdgId(),
0913              hadMothers.at(idx).status(),
0914              hadMothers.at(idx).pt(),
0915              hadMothers.at(idx).eta());
0916     }
0917   }
0918   bool hasCorrectMothers = true;
0919   if (nMothers < 1)
0920     hasCorrectMothers = false;
0921   else if (mothers.at(0) < 0)
0922     hasCorrectMothers = false;
0923   if (!hasCorrectMothers) {
0924     if (verbose)
0925       printf("    NO CORRECT MOTHER\n");
0926     return -1;
0927   }
0928   // Stopping if the particular particle has been found
0929   if (stopId >= 0 && idx == stopId)
0930     return idx;
0931
0932   // Stopping if the hadron of particular flavour has been found
0933   if (pdgId % 111111 == 0 && pdgId != 0) {
0934     if (isHadronPdgId(pdgId / 111111, hadMothers.at(idx).pdgId())) {
0935       return idx;
0936     }
0937   }
0938
0939   // Checking whether current mother satisfies selection criteria
0940   if (((hadMothers.at(idx).pdgId() == pdgId && pdgAbs == false) ||
0941        (std::abs(hadMothers.at(idx).pdgId()) == std::abs(pdgId) && pdgAbs == true)) &&
0942       (hadMothers.at(idx).status() == status || status == 0) && hasCorrectMothers) {
0943     isCorrect = true;
0944     // Adding to the list of candidates if not there and if mother of this quark has correct flavour sign
0945     const bool inList = std::find(mothChains.begin(), mothChains.end(), idx) != mothChains.end();
0946     if (!inList && mothers.at(0) >= 0 && (hadMothers.at(idx).pdgId() * pdgId > 0 || !pdgAbs)) {
0947       if (firstLast == 0 || firstLast == 1) {
0948         mothChains.push_back(idx);
0949       }
0950       if (verbose) {
0951         printf("   *");
0952       }
0953     }
0954     if (verbose) {
0955       printf("   +++");
0956     }
0957   }
0958   if (verbose) {
0959     printf("\n");
0960   }
0961
0962   if (isCorrect && firstLast == 1) {
0963     return -1;  // Stopping if only the first particle in the chain is looked for
0964   }
0965
0966   // Checking next level mothers
0967   unsigned int nDifferingMothers = 0;
0968   for (unsigned int i = 0; i < nMothers; i++) {
0969     int idx2 = mothers.at(i);
0970     if (idx2 < 0) {
0971       if (verbose)
0972         printf("^^^ Has no mother\n");
0973       continue;  // Skipping if mother's id is -1 (no mother), that means current particle is a proton
0974     }
0975     if (idx2 == idx) {
0976       if (verbose)
0977         printf("^^^ Stored as its own mother\n");
0978       continue;  // Skipping if particle is stored as its own mother
0979     }
0980     int pdg_2 = hadMothers[idx2].pdgId();
0981     // Inverting the flavour if bb oscillation detected
0982     if (isHadronPdgId(pdgId, pdg_1) && isHadronPdgId(pdgId, pdg_2) && pdg_1 * pdg_2 < 0) {
0983       pdgId *= -1;
0984       if (verbose)
0985         printf("######### Inverting flavour of the hadron\n");
0986     }
0987     // Counting how many mothers are different from this particle
0988     if ((std::abs(pdg_2) != std::abs(pdgId) && pdgAbs == true) || (pdg_2 != pdgId && pdgAbs == false)) {
0989       nDifferingMothers++;
0990     }
0991
0992     // Checking next level mother
0993     if (verbose) {
0994       printf("Checking mother %d out of %d mothers (%d -> %d), looking for pdgId: %d\n", i, nMothers, idx, idx2, pdgId);
0995     }
0996     if (firstLast == 2 && pdg_1 != pdg_2)
0997       continue;  // Requiring the same flavour when finding the last quark
0998     foundStopId = findInMothers(
0999         idx2, mothChains, hadMothersIndices, hadMothers, status, pdgId, pdgAbs, stopId, firstLast, verbose);
1000   }
1001   // Storing this particle if all its mothers are of another type and the last of its kind should be stored
1002   if (firstLast == 2 && isCorrect && nDifferingMothers >= nMothers) {
1003     if (verbose) {
1004       printf("Checking particle %d once more to store it as the last quark\n", idx);
1005     }
1006     foundStopId = findInMothers(idx, mothChains, hadMothersIndices, hadMothers, 0, pdgId, pdgAbs, stopId, 1, verbose);
1007   }
1008
1009   return foundStopId;
1010 }
1011
1012 /**
1013 * @brief Check whether a given pdgId represents neutral particle
1014 *
1015 * @param[in] pdgId
1016 * @param[in] thisParticle - a particle that is to be analysed
1017 *
1018 * @returns if the particle has a hadron among its daughters
1019 */
1020 bool GenHFHadronMatcher::isNeutralPdg(int pdgId) const {
1021   const int neutralPdgs_array[] = {9, 21, 22, 23, 25};
1022   const std::vector<int> neutralPdgs(neutralPdgs_array, neutralPdgs_array + sizeof(neutralPdgs_array) / sizeof(int));
1023   if (std::find(neutralPdgs.begin(), neutralPdgs.end(), std::abs(pdgId)) == neutralPdgs.end())
1024     return false;
1025
1026   return true;
1027 }
1028
1029 /**
1030  * Finds hadrons that have the same flavour and origin and resolve this ambiguity
1031  * @method fixExtraSameFlavours
1032  * @param  hadId                Index of the hadron being checked
1033  * @param  hadIndex             Vector of indices of each hadron pointing to the hadMothers
1034  * @param  hadMothers           Vector of gen particles (contain hadrons and all its ancestors)
1035  * @param  hadMothersIndices    Vector of indices for each particle from hadMothers
1036  * @param  isFromTopWeakDecay   Vector of values showing whether particle comes from the top weak decay
1037  * @param  LastQuarkIds         Vector of indices of last quarks for each hadron
1038  * @param  LastQuarkMotherIds   Vector of indices of mothers for each last quark from LastQuarkIds
1039  * @param  lastQuakIndices      Vector of indices pointing to particular index from the LastQuarkIds and LastQuarkMotherIds to be used for each hadron
1040  * @param  lastQuarkIndex       Index from the LastQuarkIds and LastQuarkMotherIds for this particular hadron with index hadId
1041  * @return Whether other mother with unique association has been found for the hadrons
1042  */
1043 bool GenHFHadronMatcher::fixExtraSameFlavours(const unsigned int hadId,
1044                                               const std::vector<int> &hadIndices,
1045                                               const std::vector<reco::GenParticle> &hadMothers,
1046                                               const std::vector<std::vector<int>> &hadMothersIndices,
1047                                               const std::vector<int> &isFromTopWeakDecay,
1048                                               const std::vector<std::vector<int>> &LastQuarkIds,
1049                                               const std::vector<std::vector<int>> &LastQuarkMotherIds,
1050                                               std::vector<int> &lastQuarkIndices,
1051                                               std::vector<int> &hadronFlavour,
1052                                               std::set<int> &checkedHadronIds,
1053                                               const int lastQuarkIndex) const {
1054   if (checkedHadronIds.count(hadId) != 0)
1055     return false;                  // Hadron already checked previously and should be skipped
1056   checkedHadronIds.insert(hadId);  // Putting hadron to the list of checked ones in this run
1057
1058   if (lastQuarkIndex < 0)
1059     return false;
1060   if ((int)LastQuarkIds.at(hadId).size() < lastQuarkIndex + 1)
1061     return false;
1062   int LastQuarkId = LastQuarkIds.at(hadId).at(lastQuarkIndex);
1063   int LastQuarkMotherId = LastQuarkMotherIds.at(hadId).at(lastQuarkIndex);
1064   int qmFlav = hadMothers.at(LastQuarkId).pdgId() < 0 ? -1 : 1;
1065   int hadFlavour = qmFlav * std::abs(hadMothers.at(LastQuarkMotherId).pdgId());
1066   bool ambiguityResolved = true;
1067   // If last quark has inconsistent flavour with its mother, setting the hadron flavour to gluon
1068   if ((hadMothers.at(LastQuarkId).pdgId() * hadMothers.at(LastQuarkMotherId).pdgId() < 0 &&
1069        !isNeutralPdg(hadMothers.at(LastQuarkMotherId).pdgId())) ||
1070       // If particle not coming from the Top weak decay has Top flavour
1071       (std::abs(hadronFlavour.at(hadId)) == 6 && isFromTopWeakDecay.at(hadId) == 0)) {
1072     if ((int)LastQuarkIds.at(hadId).size() > lastQuarkIndex + 1)
1073       fixExtraSameFlavours(hadId,
1074                            hadIndices,
1075                            hadMothers,
1076                            hadMothersIndices,
1077                            isFromTopWeakDecay,
1078                            LastQuarkIds,
1079                            LastQuarkMotherIds,
1080                            lastQuarkIndices,
1081                            hadronFlavour,
1082                            checkedHadronIds,
1083                            lastQuarkIndex + 1);
1084     else
1085       hadronFlavour.at(hadId) = qmFlav * 21;
1086     return true;
1087   }
1088
1089   int nSameFlavourHadrons = 0;
1090   // Looping over all previous hadrons
1091   for (unsigned int iHad = 0; iHad < hadronFlavour.size(); iHad++) {
1092     if (iHad == hadId)
1093       continue;
1094     int theLastQuarkIndex = lastQuarkIndices.at(iHad);
1095     if (theLastQuarkIndex < 0)
1096       continue;
1097     if ((int)LastQuarkMotherIds.at(iHad).size() <= theLastQuarkIndex)
1098       continue;
1099     int theLastQuarkMotherId = LastQuarkMotherIds.at(iHad).at(theLastQuarkIndex);
1100     int theHadFlavour = hadronFlavour.at(iHad);
1101     // Skipping hadrons with different flavour
1102     if (theHadFlavour == 0 || std::abs(theHadFlavour) == 21)
1103       continue;
1104     if (theHadFlavour != hadFlavour || theLastQuarkMotherId != LastQuarkMotherId)
1105       continue;
1106     ambiguityResolved = false;
1107     nSameFlavourHadrons++;
1108
1109     // Checking other b-quark mother candidates of this hadron
1110     if ((int)LastQuarkIds.at(hadId).size() > lastQuarkIndex + 1) {
1111       if (fixExtraSameFlavours(hadId,
1112                                hadIndices,
1113                                hadMothers,
1114                                hadMothersIndices,
1115                                isFromTopWeakDecay,
1116                                LastQuarkIds,
1117                                LastQuarkMotherIds,
1118                                lastQuarkIndices,
1119                                hadronFlavour,
1120                                checkedHadronIds,
1121                                lastQuarkIndex + 1)) {
1122         ambiguityResolved = true;
1123         break;
1124       }
1125     } else
1126       // Checking other b-quark mother candidates of the particular previous hadron
1127       if ((int)LastQuarkIds.at(iHad).size() > theLastQuarkIndex + 1) {
1128         if (fixExtraSameFlavours(iHad,
1129                                  hadIndices,
1130                                  hadMothers,
1131                                  hadMothersIndices,
1132                                  isFromTopWeakDecay,
1133                                  LastQuarkIds,
1134                                  LastQuarkMotherIds,
1135                                  lastQuarkIndices,
1136                                  hadronFlavour,
1137                                  checkedHadronIds,
1138                                  theLastQuarkIndex + 1)) {
1139           ambiguityResolved = true;
1140           break;
1141         };
1142       }
1143
1144   }  // End of loop over all previous hadrons
1145
1146   checkedHadronIds.erase(hadId);  // Removing the hadron from the list of checked hadrons
1147   if (nSameFlavourHadrons > 0 && !ambiguityResolved) {
1148     hadronFlavour.at(hadId) = qmFlav * 21;
1149     return true;
1150   }
1151   lastQuarkIndices.at(hadId) = lastQuarkIndex;
1152   hadronFlavour.at(hadId) = hadFlavour;
1153   return true;
1154 }
1155
1156 //define this as a plug-in
1157 DEFINE_FWK_MODULE(GenHFHadronMatcher);