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