File indexing completed on 2024-09-07 04:37:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include <memory>
0023 #include <utility>
0024 #include <vector>
0025 #include <algorithm>
0026
0027
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
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
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
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
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 }
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")},
0180 plusMothersIndicesToken_{produces<std::vector<std::vector<int>>>(
0181 "gen" + flavourStr_ + "HadPlusMothersIndices")},
0182 indexToken_{
0183 produces<std::vector<int>>("gen" + flavourStr_ + "HadIndex")},
0184 flavourToken_{produces<std::vector<int>>(
0185 "gen" + flavourStr_ +
0186 "HadFlavour")},
0187 jetIndexToken_{produces<std::vector<int>>(
0188 "gen" + flavourStr_ + "HadJetIndex")},
0189 leptonIndexToken_{produces<std::vector<int>>(
0190 "gen" + flavourStr_ +
0191 "HadLeptonIndex")},
0192 leptonHadronIndexToken_{produces<std::vector<int>>(
0193 "gen" + flavourStr_ + "HadLeptonHadronIndex")},
0194 leptonViaTauToken_{produces<std::vector<int>>(
0195 "gen" + flavourStr_ + "HadLeptonViaTau")},
0196 fromTopWeakDecayToken_{produces<std::vector<int>>(
0197 "gen" + flavourStr_ +
0198 "HadFromTopWeakDecay")},
0199 bHadronIdToken_{produces<std::vector<int>>("gen" + flavourStr_ + "HadBHadronId")}
0200
0201 {
0202
0203 }
0204
0205 GenHFHadronMatcher::~GenHFHadronMatcher() {}
0206
0207
0208
0209
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
0228
0229
0230
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
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
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
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
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
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
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
0329 if (!isHadron(flavour_, (&**hadron)))
0330 continue;
0331 if (hasHadronDaughter(flavour_, (reco::Candidate *)(&**hadron)))
0332 continue;
0333
0334 int hadronIndex = analyzeMothers((reco::Candidate *)(&**hadron),
0335 topDaughterQId,
0336 topBarDaughterQId,
0337 hadMothersCand,
0338 hadMothersIndices,
0339 nullptr,
0340 -1);
0341
0342 hadIndex.push_back(hadronIndex);
0343 hadJetIndex.push_back(jetIndex);
0344 }
0345 }
0346
0347
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
0356 if (std::find(hadMothersCand.begin(), hadMothersCand.end(), thisParticle) != hadMothersCand.end())
0357 continue;
0358
0359
0360 int hadronIndex = analyzeMothers(
0361 thisParticle, topDaughterQId, topBarDaughterQId, hadMothersCand, hadMothersIndices, nullptr, -1);
0362
0363 hadIndex.push_back(hadronIndex);
0364 hadJetIndex.push_back(-1);
0365 }
0366 }
0367
0368
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
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
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
0388 if (std::abs(leptonMother->pdgId()) == 15) {
0389 leptonViaTau = true;
0390 leptonMother = leptonMother->mother();
0391 }
0392
0393 if (leptonMother == nullptr or !isHadron(flavour_, leptonMother))
0394 continue;
0395
0396 size_t leptonHadronParticleIndex =
0397 std::find(hadMothersCand.begin(), hadMothersCand.end(), leptonMother) - hadMothersCand.begin();
0398 if (leptonHadronParticleIndex >= hadMothersCand.size())
0399 continue;
0400
0401 size_t leptonHadronIndex =
0402 std::find(hadIndex.begin(), hadIndex.end(), leptonHadronParticleIndex) - hadIndex.begin();
0403 if (leptonHadronIndex >= hadIndex.size())
0404 continue;
0405
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
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
0421 for (unsigned int hadNum = 0; hadNum < nHad; hadNum++) {
0422 unsigned int hadIdx = hadIndex.at(hadNum);
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
0431 if (hadronFlavourSign != 0) {
0432 findInMothers(
0433 hadIdx, FirstQuarkId, hadMothersIndices, hadMothers, 0, hadronFlavourSign * flavour_, false, -1, 1, false);
0434 }
0435
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
0442 for (unsigned int qId = 0; qId < FirstQuarkId.size(); qId++) {
0443
0444 const int quarkFlavourSign = flavourSign(hadMothers.at(FirstQuarkId.at(qId)).pdgId());
0445
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 }
0457
0458
0459 int hadronFlavour = 0;
0460
0461
0462 std::vector<std::pair<double, int>> lastQuark_dR_id_pairs;
0463
0464
0465 for (unsigned int qId = 0; qId < LastQuarkId.size(); qId++) {
0466 int qIdx = LastQuarkId.at(qId);
0467
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 }
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;
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);
0510
0511
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 }
0566
0567 return hadJetIndex;
0568 }
0569
0570
0571
0572
0573
0574
0575
0576
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
0597
0598
0599
0600
0601
0602
0603 bool GenHFHadronMatcher::isHadron(const int flavour, const reco::Candidate *thisParticle) const {
0604 return isHadronPdgId(flavour, thisParticle->pdgId());
0605 }
0606
0607
0608
0609
0610
0611
0612
0613
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
0624
0625
0626
0627
0628
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
0639 if (pdgId_abs / 1000 == flavour_abs)
0640 return false;
0641
0642 if (noBBbarResonances_ && pdgId_abs / 10 % 100 == 11 * flavour_abs)
0643 return false;
0644
0645 return true;
0646 }
0647
0648
0649
0650
0651
0652
0653
0654
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
0670
0671
0672
0673
0674
0675 int GenHFHadronMatcher::flavourSign(const int pdgId) const {
0676 int flavourSign = pdgId / std::abs(pdgId);
0677
0678 if (isMesonPdgId(5, pdgId))
0679 flavourSign *= -1;
0680
0681 if (pdgId % 1000 / 10 / 11 > 0)
0682 flavourSign = 0;
0683
0684 return flavourSign;
0685 }
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695 bool GenHFHadronMatcher::hasHadronDaughter(const int flavour, const reco::Candidate *thisParticle) const {
0696
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
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721
0722
0723
0724
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
0734 int hadronIndex = -1;
0735 int index = idInList(hadMothers, thisParticle);
0736 if (index < 0) {
0737 hadMothers.push_back(thisParticle);
0738 hadronIndex = hadMothers.size() - 1;
0739 } else {
0740 hadronIndex = index;
0741 }
0742
0743 int partIndex = -1;
0744 partIndex = idInList(hadMothers, thisParticle);
0745
0746
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
0760 if (isLoop) {
0761 if (prevPartIndex >= 0) {
0762 putMotherIndex(hadMothersIndices, prevPartIndex, -1);
0763 }
0764 return hadronIndex;
0765 }
0766 analyzedParticles->insert(thisParticle);
0767
0768
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;
0774 }
0775
0776
0777 if (mothIndex < 0) {
0778 hadMothers.push_back(mother);
0779 mothIndex = hadMothers.size() - 1;
0780 }
0781
0782 if (mothIndex != partIndex && partIndex >= 0) {
0783 putMotherIndex(hadMothersIndices, partIndex, mothIndex);
0784 }
0785 analyzeMothers(
0786 mother, topDaughterQId, topBarDaughterQId, hadMothers, hadMothersIndices, analyzedParticles, partIndex);
0787
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 }
0801 }
0802 }
0803
0804 analyzedParticles->erase(thisParticle);
0805
0806 if (partIndex < 0) {
0807 return hadronIndex;
0808 }
0809
0810
0811 if ((int)thisParticle->numberOfMothers() <= 0) {
0812 putMotherIndex(hadMothersIndices, partIndex, -1);
0813 }
0814
0815 return hadronIndex;
0816 }
0817
0818
0819
0820
0821
0822
0823
0824
0825
0826
0827 bool GenHFHadronMatcher::putMotherIndex(std::vector<std::vector<int>> &hadMothersIndices,
0828 int partIndex,
0829 int mothIndex) const {
0830
0831 bool inList = false;
0832 if (partIndex < 0) {
0833 return false;
0834 }
0835
0836 while ((int)hadMothersIndices.size() <= partIndex) {
0837 std::vector<int> mothersIndices;
0838 hadMothersIndices.push_back(mothersIndices);
0839 }
0840
0841 std::vector<int> *hadMotherIndices = &hadMothersIndices.at(partIndex);
0842
0843 if (mothIndex == -1) {
0844 hadMotherIndices->clear();
0845 } else {
0846
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
0856 if (!inList) {
0857 hadMotherIndices->push_back(mothIndex);
0858 }
0859
0860 return inList;
0861 }
0862
0863
0864
0865
0866
0867
0868
0869
0870
0871
0872
0873
0874
0875
0876
0877
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;
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;
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
0930 if (stopId >= 0 && idx == stopId)
0931 return idx;
0932
0933
0934 if (pdgId % 111111 == 0 && pdgId != 0) {
0935 if (isHadronPdgId(pdgId / 111111, hadMothers.at(idx).pdgId())) {
0936 return idx;
0937 }
0938 }
0939
0940
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
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;
0965 }
0966
0967
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;
0975 }
0976 if (idx2 == idx) {
0977 if (verbose)
0978 printf("^^^ Stored as its own mother\n");
0979 continue;
0980 }
0981 int pdg_2 = hadMothers[idx2].pdgId();
0982
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
0989 if ((std::abs(pdg_2) != std::abs(pdgId) && pdgAbs == true) || (pdg_2 != pdgId && pdgAbs == false)) {
0990 nDifferingMothers++;
0991 }
0992
0993
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;
0999 foundStopId = findInMothers(
1000 idx2, mothChains, hadMothersIndices, hadMothers, status, pdgId, pdgAbs, stopId, firstLast, verbose);
1001 }
1002
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
1015
1016
1017
1018
1019
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
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
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;
1057 checkedHadronIds.insert(hadId);
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
1069 if ((hadMothers.at(LastQuarkId).pdgId() * hadMothers.at(LastQuarkMotherId).pdgId() < 0 &&
1070 !isNeutralPdg(hadMothers.at(LastQuarkMotherId).pdgId())) ||
1071
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
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
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
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
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 }
1146
1147 checkedHadronIds.erase(hadId);
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
1158 DEFINE_FWK_MODULE(GenHFHadronMatcher);