Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // Collect unstable B-hadrons
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     // Collect B-hadrons, to be used in b tagging
0113     if (isBHadron(&p))
0114       bHadronIdxs.insert(-i - 1);
0115   }
0116 
0117   // Collect stable leptons and neutrinos
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;  // Skip orphans (if exists)
0129     if (p.mother()->status() == 4)
0130       continue;  // Treat particle as hadronic if directly from the incident beam (protect orphans in MINIAOD)
0131     if (isFromHadron(&p))
0132       continue;
0133     switch (absPdgId) {
0134       case 11:
0135       case 13:  // Leptons
0136       case 22:  // Photons
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   // Sort neutrinos by pT.
0148   std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
0149 
0150   // Make dressed leptons with anti-kt(0.1) algorithm
0151   //// Prepare input particle list
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   //// Run the jet algorithm
0164   fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *fjLepDef_);
0165   std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(minLeptonPt_));
0166 
0167   //// Build dressed lepton objects from the FJ output
0168   leptons->reserve(fjLepJets.size());
0169   std::set<size_t> lepDauIdxs;  // keep lepton constituents to remove from GenJet construction
0170   for (auto& fjJet : fjLepJets) {
0171     if (abs(fjJet.eta()) > maxLeptonEta_)
0172       continue;
0173 
0174     // Get jet constituents from fastJet
0175     const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
0176     // Convert to CandidatePtr
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;  // Choose one with highest pt
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     // Keep constituent indices to be used in the next step.
0206     for (auto& fjConstituent : fjConstituents) {
0207       lepDauIdxs.insert(fjConstituent.user_index());
0208     }
0209   }
0210 
0211   // Now proceed to jets.
0212   // Jets: anti-kt excluding the e, mu, nu, and photons in selected leptons.
0213   //// Prepare input particle list. Remove particles used in lepton clusters, neutrinos
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   //// Also don't forget to put B hadrons
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   //// Run the jet algorithm
0244   fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *fjJetDef_);
0245   std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(minJetPt_));
0246 
0247   /// Build jets
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     // Get jet constituents from fastJet
0255     const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
0256     // Convert to CandidatePtr
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   // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
0286   // NOTE : A C++ trick, use do-while instead of long-nested if-statements.
0287   do {
0288     if (bjetIdxs.size() < 2)
0289       break;
0290 
0291     // Note : we will do dilepton or semilepton channel only
0292     if (leptons->size() == 2 and neutrinos->size() >= 2) {
0293       // Start from dilepton channel
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) {  // Consider leading 2 neutrinos. Neutrino vector is already sorted by pT
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       // Contiue to top quarks
0335       dm = 1e9;  // Reset once again for top combination.
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       // Recalculate lepton charges (needed due to initialization of leptons wrt sign of a charge)
0361       const int lepQ1 = lepton1.charge();
0362       const int lepQ2 = lepton2.charge();
0363 
0364       // Put all of them into candidate collection
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       // Then continue to the semi-leptonic channel
0391       const auto& lepton = leptons->at(0);
0392       if (lepton.pt() < minLeptonPtSemilepton_ or std::abs(lepton.eta()) > maxLeptonEtaSemilepton_)
0393         break;
0394 
0395       // Skip event if there are veto leptons
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       // Calculate MET
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       // Solve pz to give wMass_^2 = (lepE+energy)^2 - (lepPx+metX)^2 - (lepPy+metY)^2 - (lepPz+pz)^2
0421       // -> (wMass_^2)/2 = lepE*sqrt(metPt^2+pz^2) - lepPx*metX - lepPy*metY - lepPz*pz
0422       // -> C := (wMass_^2)/2 + lepPx*metX + lepPy*metY
0423       // -> lepE^2*(metPt^2+pz^2) = C^2 + 2*C*lepPz*pz + lepPz^2*pz^2
0424       // -> (lepE^2-lepPz^2)*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
0425       // -> lepPt^2*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
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       // Continue to build leptonic pseudo top
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       // Build hadronic pseudo W, take leading two jets (ljetIdxs are (implicitly) sorted by pT)
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       // Contiue to top quarks
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       // Put all of them into candidate collection
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)  // If pseudtop decay tree is completed
0501   {
0502     // t->W+b
0503     pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 2));  // t->W
0504     pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 3));  // t->b
0505     pseudoTop->at(2).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0));    // t->W
0506     pseudoTop->at(3).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0));    // t->b
0507 
0508     // W->lv or W->jj
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     // tbar->W-b
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     // W->jj
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;  // Skip incident beam
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   // Do not consider this particle if it has B hadron daughter
0566   // For example, B* -> B0 + photon; then we drop B* and take B0 only
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;  // Fundamental particles and MC internals
0579   if (absPdgId >= 1000000000)
0580     return false;  // Nuclei, +-10LZZZAAAI
0581 
0582   // General form of PDG ID is 7 digit form
0583   // +- n nr nL nq1 nq2 nq3 nJ
0584   //const int nJ = absPdgId % 10; // Spin
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;  // Diquarks
0591   if (nq1 == 0 and nq2 == 5)
0592     return true;  // B mesons
0593   if (nq1 == 5)
0594     return true;  // B baryons
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!