File indexing completed on 2022-01-25 05:52:08
0001 #include "DataFormats/Math/interface/deltaR.h"
0002 #include "CommonTools/Utils/interface/PtComparator.h"
0003
0004 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0005 #include "fastjet/ClusterSequence.hh"
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 #include "DataFormats/Candidate/interface/Candidate.h"
0013 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0014
0015 #include "fastjet/JetDefinition.hh"
0016 #include <set>
0017
0018 class PseudoTopProducer : public edm::stream::EDProducer<> {
0019 public:
0020 PseudoTopProducer(const edm::ParameterSet& pset);
0021 void produce(edm::Event& event, const edm::EventSetup& eventSetup) override;
0022
0023 private:
0024 bool isFromHadron(const reco::Candidate* p) const;
0025 bool isBHadron(const reco::Candidate* p) const;
0026 bool isBHadron(const unsigned int pdgId) const;
0027 void insertAllDaughters(const reco::Candidate* p, std::set<const reco::Candidate*>& list) const;
0028
0029 const reco::Candidate* getLast(const reco::Candidate* p);
0030 reco::GenParticleRef buildGenParticle(const reco::Candidate* p,
0031 reco::GenParticleRefProd& refHandle,
0032 std::auto_ptr<reco::GenParticleCollection>& outColl) const;
0033 typedef reco::Particle::LorentzVector LorentzVector;
0034
0035 private:
0036 const edm::EDGetTokenT<edm::View<reco::Candidate> > finalStateToken_;
0037 const edm::EDGetTokenT<edm::View<reco::Candidate> > genParticleToken_;
0038 const double minLeptonPt_, maxLeptonEta_, minJetPt_, maxJetEta_;
0039 const double wMass_, tMass_;
0040 const double minLeptonPtDilepton_, maxLeptonEtaDilepton_;
0041 const double minDileptonMassDilepton_;
0042 const double minLeptonPtSemilepton_, maxLeptonEtaSemilepton_;
0043 const double minVetoLeptonPtSemilepton_, maxVetoLeptonEtaSemilepton_;
0044 const double minMETSemiLepton_, minMtWSemiLepton_;
0045
0046 typedef fastjet::JetDefinition JetDef;
0047 std::shared_ptr<JetDef> fjLepDef_, fjJetDef_;
0048 reco::Particle::Point genVertex_;
0049 };
0050
0051 using namespace std;
0052 using namespace edm;
0053 using namespace reco;
0054
0055 PseudoTopProducer::PseudoTopProducer(const edm::ParameterSet& pset)
0056 : finalStateToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("finalStates"))),
0057 genParticleToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("genParticles"))),
0058 minLeptonPt_(pset.getParameter<double>("minLeptonPt")),
0059 maxLeptonEta_(pset.getParameter<double>("maxLeptonEta")),
0060 minJetPt_(pset.getParameter<double>("minJetPt")),
0061 maxJetEta_(pset.getParameter<double>("maxJetEta")),
0062 wMass_(pset.getParameter<double>("wMass")),
0063 tMass_(pset.getParameter<double>("tMass")),
0064 minLeptonPtDilepton_(pset.getParameter<double>("minLeptonPtDilepton")),
0065 maxLeptonEtaDilepton_(pset.getParameter<double>("maxLeptonEtaDilepton")),
0066 minDileptonMassDilepton_(pset.getParameter<double>("minDileptonMassDilepton")),
0067 minLeptonPtSemilepton_(pset.getParameter<double>("minLeptonPtSemilepton")),
0068 maxLeptonEtaSemilepton_(pset.getParameter<double>("maxLeptonEtaSemilepton")),
0069 minVetoLeptonPtSemilepton_(pset.getParameter<double>("minVetoLeptonPtSemilepton")),
0070 maxVetoLeptonEtaSemilepton_(pset.getParameter<double>("maxVetoLeptonEtaSemilepton")),
0071 minMETSemiLepton_(pset.getParameter<double>("minMETSemiLepton")),
0072 minMtWSemiLepton_(pset.getParameter<double>("minMtWSemiLepton")) {
0073 const double leptonConeSize = pset.getParameter<double>("leptonConeSize");
0074 const double jetConeSize = pset.getParameter<double>("jetConeSize");
0075 fjLepDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, leptonConeSize);
0076 fjJetDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, jetConeSize);
0077
0078 genVertex_ = reco::Particle::Point(0, 0, 0);
0079
0080 produces<reco::GenParticleCollection>("neutrinos");
0081 produces<reco::GenJetCollection>("leptons");
0082 produces<reco::GenJetCollection>("jets");
0083
0084 produces<reco::GenParticleCollection>();
0085 }
0086
0087 void PseudoTopProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup) {
0088 edm::Handle<edm::View<reco::Candidate> > finalStateHandle;
0089 event.getByToken(finalStateToken_, finalStateHandle);
0090
0091 edm::Handle<edm::View<reco::Candidate> > genParticleHandle;
0092 event.getByToken(genParticleToken_, genParticleHandle);
0093
0094 std::unique_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
0095 std::unique_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
0096 std::unique_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
0097 auto neutrinosRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("neutrinos");
0098 auto leptonsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("leptons");
0099 auto jetsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("jets");
0100
0101 std::unique_ptr<reco::GenParticleCollection> pseudoTop(new reco::GenParticleCollection);
0102 auto pseudoTopRefHandle = event.getRefBeforePut<reco::GenParticleCollection>();
0103
0104
0105 std::set<int> bHadronIdxs;
0106 for (size_t i = 0, n = genParticleHandle->size(); i < n; ++i) {
0107 const reco::Candidate& p = genParticleHandle->at(i);
0108 const int status = p.status();
0109 if (status == 1)
0110 continue;
0111
0112
0113 if (isBHadron(&p))
0114 bHadronIdxs.insert(-i - 1);
0115 }
0116
0117
0118 size_t nStables = 0;
0119 std::vector<size_t> leptonIdxs;
0120 for (size_t i = 0, n = finalStateHandle->size(); i < n; ++i) {
0121 const reco::Candidate& p = finalStateHandle->at(i);
0122 const int absPdgId = abs(p.pdgId());
0123 if (p.status() != 1)
0124 continue;
0125
0126 ++nStables;
0127 if (p.numberOfMothers() == 0)
0128 continue;
0129 if (p.mother()->status() == 4)
0130 continue;
0131 if (isFromHadron(&p))
0132 continue;
0133 switch (absPdgId) {
0134 case 11:
0135 case 13:
0136 case 22:
0137 leptonIdxs.push_back(i);
0138 break;
0139 case 12:
0140 case 14:
0141 case 16:
0142 neutrinos->push_back(reco::GenParticle(p.charge(), p.p4(), p.vertex(), p.pdgId(), p.status(), true));
0143 break;
0144 }
0145 }
0146
0147
0148 std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
0149
0150
0151
0152 std::vector<fastjet::PseudoJet> fjLepInputs;
0153 fjLepInputs.reserve(leptonIdxs.size());
0154 for (auto index : leptonIdxs) {
0155 const reco::Candidate& p = finalStateHandle->at(index);
0156 if (std::isnan(p.pt()) or p.pt() <= 0)
0157 continue;
0158
0159 fjLepInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
0160 fjLepInputs.back().set_user_index(index);
0161 }
0162
0163
0164 fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *fjLepDef_);
0165 std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(minLeptonPt_));
0166
0167
0168 leptons->reserve(fjLepJets.size());
0169 std::set<size_t> lepDauIdxs;
0170 for (auto& fjJet : fjLepJets) {
0171 if (abs(fjJet.eta()) > maxLeptonEta_)
0172 continue;
0173
0174
0175 const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
0176
0177 std::vector<reco::CandidatePtr> constituents;
0178 reco::CandidatePtr lepCand;
0179 for (auto& fjConstituent : fjConstituents) {
0180 const size_t index = fjConstituent.user_index();
0181 reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
0182 const int absPdgId = abs(cand->pdgId());
0183 if (absPdgId == 11 or absPdgId == 13) {
0184 if (lepCand.isNonnull() and lepCand->pt() > cand->pt())
0185 continue;
0186 lepCand = cand;
0187 }
0188 constituents.push_back(cand);
0189 }
0190 if (lepCand.isNull())
0191 continue;
0192
0193 const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
0194 reco::GenJet lepJet;
0195 reco::writeSpecific(lepJet, jetP4, genVertex_, constituents);
0196
0197 lepJet.setPdgId(lepCand->pdgId());
0198 lepJet.setCharge(lepCand->charge());
0199
0200 const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
0201 lepJet.setJetArea(jetArea);
0202
0203 leptons->push_back(lepJet);
0204
0205
0206 for (auto& fjConstituent : fjConstituents) {
0207 lepDauIdxs.insert(fjConstituent.user_index());
0208 }
0209 }
0210
0211
0212
0213
0214 std::vector<fastjet::PseudoJet> fjJetInputs;
0215 fjJetInputs.reserve(nStables);
0216 for (size_t i = 0, n = finalStateHandle->size(); i < n; ++i) {
0217 const reco::Candidate& p = finalStateHandle->at(i);
0218 if (p.status() != 1)
0219 continue;
0220 if (std::isnan(p.pt()) or p.pt() <= 0)
0221 continue;
0222
0223 const int absId = std::abs(p.pdgId());
0224 if (absId == 12 or absId == 14 or absId == 16)
0225 continue;
0226 if (lepDauIdxs.find(i) != lepDauIdxs.end())
0227 continue;
0228
0229 fjJetInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
0230 fjJetInputs.back().set_user_index(i);
0231 }
0232
0233 for (auto index : bHadronIdxs) {
0234 const reco::Candidate& p = genParticleHandle->at(abs(index + 1));
0235 if (std::isnan(p.pt()) or p.pt() <= 0)
0236 continue;
0237
0238 const double scale = 1e-20 / p.p();
0239 fjJetInputs.push_back(fastjet::PseudoJet(p.px() * scale, p.py() * scale, p.pz() * scale, p.energy() * scale));
0240 fjJetInputs.back().set_user_index(index);
0241 }
0242
0243
0244 fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *fjJetDef_);
0245 std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(minJetPt_));
0246
0247
0248 jets->reserve(fjJets.size());
0249 std::vector<size_t> bjetIdxs, ljetIdxs;
0250 for (auto& fjJet : fjJets) {
0251 if (abs(fjJet.eta()) > maxJetEta_)
0252 continue;
0253
0254
0255 const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
0256
0257 std::vector<reco::CandidatePtr> constituents;
0258 bool hasBHadron = false;
0259 for (size_t j = 0, m = fjConstituents.size(); j < m; ++j) {
0260 const int index = fjConstituents[j].user_index();
0261 if (bHadronIdxs.find(index) != bHadronIdxs.end()) {
0262 hasBHadron = true;
0263 continue;
0264 }
0265 reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
0266 constituents.push_back(cand);
0267 }
0268
0269 const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
0270 reco::GenJet genJet;
0271 reco::writeSpecific(genJet, jetP4, genVertex_, constituents);
0272
0273 const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
0274 genJet.setJetArea(jetArea);
0275 if (hasBHadron) {
0276 genJet.setPdgId(5);
0277 bjetIdxs.push_back(jets->size());
0278 } else {
0279 ljetIdxs.push_back(jets->size());
0280 }
0281
0282 jets->push_back(genJet);
0283 }
0284
0285
0286
0287 do {
0288 if (bjetIdxs.size() < 2)
0289 break;
0290
0291
0292 if (leptons->size() == 2 and neutrinos->size() >= 2) {
0293
0294 const int q1 = leptons->at(0).charge();
0295 const int q2 = leptons->at(1).charge();
0296 if (q1 * q2 > 0)
0297 break;
0298
0299 const auto& lepton1 = q1 > 0 ? leptons->at(0) : leptons->at(1);
0300 const auto& lepton2 = q1 > 0 ? leptons->at(1) : leptons->at(0);
0301
0302 if (lepton1.pt() < minLeptonPtDilepton_ or std::abs(lepton1.eta()) > maxLeptonEtaDilepton_)
0303 break;
0304 if (lepton2.pt() < minLeptonPtDilepton_ or std::abs(lepton2.eta()) > maxLeptonEtaDilepton_)
0305 break;
0306 if ((lepton1.p4() + lepton2.p4()).mass() < minDileptonMassDilepton_)
0307 break;
0308
0309 double dm = 1e9;
0310 int selNu1 = -1, selNu2 = -1;
0311 for (int i = 0; i < 2; ++i) {
0312 const double dm1 = std::abs((lepton1.p4() + neutrinos->at(i).p4()).mass() - wMass_);
0313 for (int j = 0; j < 2; ++j) {
0314 if (i == j)
0315 continue;
0316 const double dm2 = std::abs((lepton2.p4() + neutrinos->at(j).p4()).mass() - wMass_);
0317 const double newDm = dm1 + dm2;
0318
0319 if (newDm < dm) {
0320 dm = newDm;
0321 selNu1 = i;
0322 selNu2 = j;
0323 }
0324 }
0325 }
0326 if (dm >= 1e9)
0327 break;
0328
0329 const auto& nu1 = neutrinos->at(selNu1);
0330 const auto& nu2 = neutrinos->at(selNu2);
0331 const auto w1LVec = lepton1.p4() + nu1.p4();
0332 const auto w2LVec = lepton2.p4() + nu2.p4();
0333
0334
0335 dm = 1e9;
0336 int selB1 = -1, selB2 = -1;
0337 for (auto i : bjetIdxs) {
0338 const double dm1 = std::abs((w1LVec + jets->at(i).p4()).mass() - tMass_);
0339 for (auto j : bjetIdxs) {
0340 if (i == j)
0341 continue;
0342 const double dm2 = std::abs((w2LVec + jets->at(j).p4()).mass() - tMass_);
0343 const double newDm = dm1 + dm2;
0344
0345 if (newDm < dm) {
0346 dm = newDm;
0347 selB1 = i;
0348 selB2 = j;
0349 }
0350 }
0351 }
0352 if (dm >= 1e9)
0353 break;
0354
0355 const auto& bJet1 = jets->at(selB1);
0356 const auto& bJet2 = jets->at(selB2);
0357 const auto t1LVec = w1LVec + bJet1.p4();
0358 const auto t2LVec = w2LVec + bJet2.p4();
0359
0360
0361 const int lepQ1 = lepton1.charge();
0362 const int lepQ2 = lepton2.charge();
0363
0364
0365 reco::GenParticle t1(lepQ1 * 2 / 3., t1LVec, genVertex_, lepQ1 * 6, 3, false);
0366 reco::GenParticle w1(lepQ1, w1LVec, genVertex_, lepQ1 * 24, 3, true);
0367 reco::GenParticle b1(-lepQ1 / 3., bJet1.p4(), genVertex_, lepQ1 * 5, 1, true);
0368 reco::GenParticle l1(lepQ1, lepton1.p4(), genVertex_, lepton1.pdgId(), 1, true);
0369 reco::GenParticle n1(0, nu1.p4(), genVertex_, nu1.pdgId(), 1, true);
0370
0371 reco::GenParticle t2(lepQ2 * 2 / 3., t2LVec, genVertex_, lepQ2 * 6, 3, false);
0372 reco::GenParticle w2(lepQ2, w2LVec, genVertex_, lepQ2 * 24, 3, true);
0373 reco::GenParticle b2(-lepQ2 / 3., bJet2.p4(), genVertex_, lepQ2 * 5, 1, true);
0374 reco::GenParticle l2(lepQ2, lepton2.p4(), genVertex_, lepton2.pdgId(), 1, true);
0375 reco::GenParticle n2(0, nu2.p4(), genVertex_, nu2.pdgId(), 1, true);
0376
0377 pseudoTop->push_back(t1);
0378 pseudoTop->push_back(t2);
0379
0380 pseudoTop->push_back(w1);
0381 pseudoTop->push_back(b1);
0382 pseudoTop->push_back(l1);
0383 pseudoTop->push_back(n1);
0384
0385 pseudoTop->push_back(w2);
0386 pseudoTop->push_back(b2);
0387 pseudoTop->push_back(l2);
0388 pseudoTop->push_back(n2);
0389 } else if (!leptons->empty() and !neutrinos->empty() and ljetIdxs.size() >= 2) {
0390
0391 const auto& lepton = leptons->at(0);
0392 if (lepton.pt() < minLeptonPtSemilepton_ or std::abs(lepton.eta()) > maxLeptonEtaSemilepton_)
0393 break;
0394
0395
0396 bool hasVetoLepton = false;
0397 for (auto vetoLeptonItr = next(leptons->begin()); vetoLeptonItr != leptons->end(); ++vetoLeptonItr) {
0398 if (vetoLeptonItr->pt() > minVetoLeptonPtSemilepton_ and
0399 std::abs(vetoLeptonItr->eta()) < maxVetoLeptonEtaSemilepton_) {
0400 hasVetoLepton = true;
0401 }
0402 }
0403 if (hasVetoLepton)
0404 break;
0405
0406
0407 double metX = 0, metY = 0;
0408 for (const auto& neutrino : *neutrinos) {
0409 metX += neutrino.px();
0410 metY += neutrino.py();
0411 }
0412 const double metPt = std::hypot(metX, metY);
0413 if (metPt < minMETSemiLepton_)
0414 break;
0415
0416 const double mtW = std::sqrt(2 * lepton.pt() * metPt * cos(lepton.phi() - atan2(metX, metY)));
0417 if (mtW < minMtWSemiLepton_)
0418 break;
0419
0420
0421
0422
0423
0424
0425
0426 const double lpz = lepton.pz(), le = lepton.energy(), lpt = lepton.pt();
0427 const double cf = (wMass_ * wMass_) / 2 + lepton.px() * metX + lepton.py() * metY;
0428 const double det = cf * cf * lpz * lpz - lpt * lpt * (le * le * metPt * metPt - cf * cf);
0429 const double pz = (cf * lpz + (cf < 0 ? -sqrt(det) : sqrt(det))) / lpt / lpt;
0430 const reco::Candidate::LorentzVector nu1P4(metX, metY, pz, std::hypot(metPt, pz));
0431 const auto w1LVec = lepton.p4() + nu1P4;
0432
0433
0434 double minDR = 1e9;
0435 int selB1 = -1;
0436 for (auto b1Itr = bjetIdxs.begin(); b1Itr != bjetIdxs.end(); ++b1Itr) {
0437 const double dR = deltaR(jets->at(*b1Itr).p4(), w1LVec);
0438 if (dR < minDR) {
0439 selB1 = *b1Itr;
0440 minDR = dR;
0441 }
0442 }
0443 if (selB1 == -1)
0444 break;
0445 const auto& bJet1 = jets->at(selB1);
0446 const auto t1LVec = w1LVec + bJet1.p4();
0447
0448
0449 const auto& wJet1 = jets->at(ljetIdxs[0]);
0450 const auto& wJet2 = jets->at(ljetIdxs[1]);
0451 const auto w2LVec = wJet1.p4() + wJet2.p4();
0452
0453
0454 double dm = 1e9;
0455 int selB2 = -1;
0456 for (auto i : bjetIdxs) {
0457 if (int(i) == selB1)
0458 continue;
0459 const double newDm = std::abs((w2LVec + jets->at(i).p4()).mass() - tMass_);
0460 if (newDm < dm) {
0461 dm = newDm;
0462 selB2 = i;
0463 }
0464 }
0465 if (dm >= 1e9)
0466 break;
0467
0468 const auto& bJet2 = jets->at(selB2);
0469 const auto t2LVec = w2LVec + bJet2.p4();
0470
0471 const int q = lepton.charge();
0472
0473 reco::GenParticle t1(q * 2 / 3., t1LVec, genVertex_, q * 6, 3, false);
0474 reco::GenParticle w1(q, w1LVec, genVertex_, q * 24, 3, true);
0475 reco::GenParticle b1(-q / 3., bJet1.p4(), genVertex_, q * 5, 1, true);
0476 reco::GenParticle l1(q, lepton.p4(), genVertex_, lepton.pdgId(), 1, true);
0477 reco::GenParticle n1(0, nu1P4, genVertex_, q * (std::abs(lepton.pdgId()) + 1), 1, true);
0478
0479 reco::GenParticle t2(-q * 2 / 3., t2LVec, genVertex_, -q * 6, 3, false);
0480 reco::GenParticle w2(-q, w2LVec, genVertex_, -q * 24, 3, true);
0481 reco::GenParticle b2(0, bJet2.p4(), genVertex_, -q * 5, 1, true);
0482 reco::GenParticle u2(0, wJet1.p4(), genVertex_, -2 * q, 1, true);
0483 reco::GenParticle d2(0, wJet2.p4(), genVertex_, q, 1, true);
0484
0485 pseudoTop->push_back(t1);
0486 pseudoTop->push_back(t2);
0487
0488 pseudoTop->push_back(w1);
0489 pseudoTop->push_back(b1);
0490 pseudoTop->push_back(l1);
0491 pseudoTop->push_back(n1);
0492
0493 pseudoTop->push_back(w2);
0494 pseudoTop->push_back(b2);
0495 pseudoTop->push_back(u2);
0496 pseudoTop->push_back(d2);
0497 }
0498 } while (false);
0499
0500 if (pseudoTop->size() == 10)
0501 {
0502
0503 pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 2));
0504 pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 3));
0505 pseudoTop->at(2).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0));
0506 pseudoTop->at(3).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0));
0507
0508
0509 pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 4));
0510 pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 5));
0511 pseudoTop->at(4).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
0512 pseudoTop->at(5).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
0513
0514
0515 pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 6));
0516 pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 7));
0517 pseudoTop->at(6).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
0518 pseudoTop->at(7).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
0519
0520
0521 pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 8));
0522 pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 9));
0523 pseudoTop->at(8).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
0524 pseudoTop->at(9).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
0525 }
0526
0527 event.put(std::move(neutrinos), "neutrinos");
0528 event.put(std::move(leptons), "leptons");
0529 event.put(std::move(jets), "jets");
0530
0531 event.put(std::move(pseudoTop));
0532 }
0533
0534 const reco::Candidate* PseudoTopProducer::getLast(const reco::Candidate* p) {
0535 for (size_t i = 0, n = p->numberOfDaughters(); i < n; ++i) {
0536 const reco::Candidate* dau = p->daughter(i);
0537 if (p->pdgId() == dau->pdgId())
0538 return getLast(dau);
0539 }
0540 return p;
0541 }
0542
0543 bool PseudoTopProducer::isFromHadron(const reco::Candidate* p) const {
0544 for (size_t i = 0, n = p->numberOfMothers(); i < n; ++i) {
0545 const reco::Candidate* mother = p->mother(i);
0546 if (!mother)
0547 continue;
0548 if (mother->status() == 4 or mother->numberOfMothers() == 0)
0549 continue;
0550 const int pdgId = abs(mother->pdgId());
0551
0552 if (pdgId > 100)
0553 return true;
0554 else if (isFromHadron(mother))
0555 return true;
0556 }
0557 return false;
0558 }
0559
0560 bool PseudoTopProducer::isBHadron(const reco::Candidate* p) const {
0561 const unsigned int absPdgId = abs(p->pdgId());
0562 if (!isBHadron(absPdgId))
0563 return false;
0564
0565
0566
0567 for (int i = 0, n = p->numberOfDaughters(); i < n; ++i) {
0568 const reco::Candidate* dau = p->daughter(i);
0569 if (isBHadron(abs(dau->pdgId())))
0570 return false;
0571 }
0572
0573 return true;
0574 }
0575
0576 bool PseudoTopProducer::isBHadron(const unsigned int absPdgId) const {
0577 if (absPdgId <= 100)
0578 return false;
0579 if (absPdgId >= 1000000000)
0580 return false;
0581
0582
0583
0584
0585 const int nq3 = (absPdgId / 10) % 10;
0586 const int nq2 = (absPdgId / 100) % 10;
0587 const int nq1 = (absPdgId / 1000) % 10;
0588
0589 if (nq3 == 0)
0590 return false;
0591 if (nq1 == 0 and nq2 == 5)
0592 return true;
0593 if (nq1 == 5)
0594 return true;
0595
0596 return false;
0597 }
0598
0599 reco::GenParticleRef PseudoTopProducer::buildGenParticle(const reco::Candidate* p,
0600 reco::GenParticleRefProd& refHandle,
0601 std::auto_ptr<reco::GenParticleCollection>& outColl) const {
0602 reco::GenParticle pOut(*dynamic_cast<const reco::GenParticle*>(p));
0603 pOut.clearMothers();
0604 pOut.clearDaughters();
0605 pOut.resetMothers(refHandle.id());
0606 pOut.resetDaughters(refHandle.id());
0607
0608 outColl->push_back(pOut);
0609
0610 return reco::GenParticleRef(refHandle, outColl->size() - 1);
0611 }
0612
0613 #include "FWCore/Framework/interface/MakerMacros.h"
0614 DEFINE_FWK_MODULE(PseudoTopProducer);