Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:13 

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 
0028   typedef reco::Particle::LorentzVector LorentzVector;
0029 
0030 private:
0031   const edm::EDGetTokenT<edm::View<reco::Candidate> > finalStateToken_;
0032   const edm::EDGetTokenT<edm::View<reco::Candidate> > genParticleToken_;
0033   const double minLeptonPt_, maxLeptonEta_, minJetPt_, maxJetEta_;
0034   const double wMass_, tMass_;
0035   const double minLeptonPtDilepton_, maxLeptonEtaDilepton_;
0036   const double minDileptonMassDilepton_;
0037   const double minLeptonPtSemilepton_, maxLeptonEtaSemilepton_;
0038   const double minVetoLeptonPtSemilepton_, maxVetoLeptonEtaSemilepton_;
0039   const double minMETSemiLepton_, minMtWSemiLepton_;
0040 
0041   typedef fastjet::JetDefinition JetDef;
0042   std::shared_ptr<JetDef> fjLepDef_, fjJetDef_;
0043   reco::Particle::Point genVertex_;
0044 };
0045 
0046 using namespace std;
0047 using namespace edm;
0048 using namespace reco;
0049 
0050 PseudoTopProducer::PseudoTopProducer(const edm::ParameterSet& pset)
0051     : finalStateToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("finalStates"))),
0052       genParticleToken_(consumes<edm::View<reco::Candidate> >(pset.getParameter<edm::InputTag>("genParticles"))),
0053       minLeptonPt_(pset.getParameter<double>("minLeptonPt")),
0054       maxLeptonEta_(pset.getParameter<double>("maxLeptonEta")),
0055       minJetPt_(pset.getParameter<double>("minJetPt")),
0056       maxJetEta_(pset.getParameter<double>("maxJetEta")),
0057       wMass_(pset.getParameter<double>("wMass")),
0058       tMass_(pset.getParameter<double>("tMass")),
0059       minLeptonPtDilepton_(pset.getParameter<double>("minLeptonPtDilepton")),
0060       maxLeptonEtaDilepton_(pset.getParameter<double>("maxLeptonEtaDilepton")),
0061       minDileptonMassDilepton_(pset.getParameter<double>("minDileptonMassDilepton")),
0062       minLeptonPtSemilepton_(pset.getParameter<double>("minLeptonPtSemilepton")),
0063       maxLeptonEtaSemilepton_(pset.getParameter<double>("maxLeptonEtaSemilepton")),
0064       minVetoLeptonPtSemilepton_(pset.getParameter<double>("minVetoLeptonPtSemilepton")),
0065       maxVetoLeptonEtaSemilepton_(pset.getParameter<double>("maxVetoLeptonEtaSemilepton")),
0066       minMETSemiLepton_(pset.getParameter<double>("minMETSemiLepton")),
0067       minMtWSemiLepton_(pset.getParameter<double>("minMtWSemiLepton")) {
0068   const double leptonConeSize = pset.getParameter<double>("leptonConeSize");
0069   const double jetConeSize = pset.getParameter<double>("jetConeSize");
0070   fjLepDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, leptonConeSize);
0071   fjJetDef_ = std::make_shared<JetDef>(fastjet::antikt_algorithm, jetConeSize);
0072 
0073   genVertex_ = reco::Particle::Point(0, 0, 0);
0074 
0075   produces<reco::GenParticleCollection>("neutrinos");
0076   produces<reco::GenJetCollection>("leptons");
0077   produces<reco::GenJetCollection>("jets");
0078 
0079   produces<reco::GenParticleCollection>();
0080 }
0081 
0082 void PseudoTopProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup) {
0083   edm::Handle<edm::View<reco::Candidate> > finalStateHandle;
0084   event.getByToken(finalStateToken_, finalStateHandle);
0085 
0086   edm::Handle<edm::View<reco::Candidate> > genParticleHandle;
0087   event.getByToken(genParticleToken_, genParticleHandle);
0088 
0089   std::unique_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
0090   std::unique_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
0091   std::unique_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
0092   auto neutrinosRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("neutrinos");
0093   auto leptonsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("leptons");
0094   auto jetsRefHandle = event.getRefBeforePut<reco::GenJetCollection>("jets");
0095 
0096   std::unique_ptr<reco::GenParticleCollection> pseudoTop(new reco::GenParticleCollection);
0097   auto pseudoTopRefHandle = event.getRefBeforePut<reco::GenParticleCollection>();
0098 
0099   // Collect unstable B-hadrons
0100   std::set<int> bHadronIdxs;
0101   for (size_t i = 0, n = genParticleHandle->size(); i < n; ++i) {
0102     const reco::Candidate& p = genParticleHandle->at(i);
0103     const int status = p.status();
0104     if (status == 1)
0105       continue;
0106 
0107     // Collect B-hadrons, to be used in b tagging
0108     if (isBHadron(&p))
0109       bHadronIdxs.insert(-i - 1);
0110   }
0111 
0112   // Collect stable leptons and neutrinos
0113   size_t nStables = 0;
0114   std::vector<size_t> leptonIdxs;
0115   for (size_t i = 0, n = finalStateHandle->size(); i < n; ++i) {
0116     const reco::Candidate& p = finalStateHandle->at(i);
0117     const int absPdgId = abs(p.pdgId());
0118     if (p.status() != 1)
0119       continue;
0120 
0121     ++nStables;
0122     if (p.numberOfMothers() == 0)
0123       continue;  // Skip orphans (if exists)
0124     if (p.mother()->status() == 4)
0125       continue;  // Treat particle as hadronic if directly from the incident beam (protect orphans in MINIAOD)
0126     if (isFromHadron(&p))
0127       continue;
0128     switch (absPdgId) {
0129       case 11:
0130       case 13:  // Leptons
0131       case 22:  // Photons
0132         leptonIdxs.push_back(i);
0133         break;
0134       case 12:
0135       case 14:
0136       case 16:
0137         neutrinos->push_back(reco::GenParticle(p.charge(), p.p4(), p.vertex(), p.pdgId(), p.status(), true));
0138         break;
0139     }
0140   }
0141 
0142   // Sort neutrinos by pT.
0143   std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
0144 
0145   // Make dressed leptons with anti-kt(0.1) algorithm
0146   //// Prepare input particle list
0147   std::vector<fastjet::PseudoJet> fjLepInputs;
0148   fjLepInputs.reserve(leptonIdxs.size());
0149   for (auto index : leptonIdxs) {
0150     const reco::Candidate& p = finalStateHandle->at(index);
0151     if (std::isnan(p.pt()) or p.pt() <= 0)
0152       continue;
0153 
0154     fjLepInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
0155     fjLepInputs.back().set_user_index(index);
0156   }
0157 
0158   //// Run the jet algorithm
0159   fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *fjLepDef_);
0160   std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(minLeptonPt_));
0161 
0162   //// Build dressed lepton objects from the FJ output
0163   leptons->reserve(fjLepJets.size());
0164   std::set<size_t> lepDauIdxs;  // keep lepton constituents to remove from GenJet construction
0165   for (auto& fjJet : fjLepJets) {
0166     if (abs(fjJet.eta()) > maxLeptonEta_)
0167       continue;
0168 
0169     // Get jet constituents from fastJet
0170     const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
0171     // Convert to CandidatePtr
0172     std::vector<reco::CandidatePtr> constituents;
0173     reco::CandidatePtr lepCand;
0174     for (auto& fjConstituent : fjConstituents) {
0175       const size_t index = fjConstituent.user_index();
0176       reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
0177       const int absPdgId = abs(cand->pdgId());
0178       if (absPdgId == 11 or absPdgId == 13) {
0179         if (lepCand.isNonnull() and lepCand->pt() > cand->pt())
0180           continue;  // Choose one with highest pt
0181         lepCand = cand;
0182       }
0183       constituents.push_back(cand);
0184     }
0185     if (lepCand.isNull())
0186       continue;
0187 
0188     const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
0189     reco::GenJet lepJet;
0190     reco::writeSpecific(lepJet, jetP4, genVertex_, constituents);
0191 
0192     lepJet.setPdgId(lepCand->pdgId());
0193     lepJet.setCharge(lepCand->charge());
0194 
0195     const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
0196     lepJet.setJetArea(jetArea);
0197 
0198     leptons->push_back(lepJet);
0199 
0200     // Keep constituent indices to be used in the next step.
0201     for (auto& fjConstituent : fjConstituents) {
0202       lepDauIdxs.insert(fjConstituent.user_index());
0203     }
0204   }
0205 
0206   // Now proceed to jets.
0207   // Jets: anti-kt excluding the e, mu, nu, and photons in selected leptons.
0208   //// Prepare input particle list. Remove particles used in lepton clusters, neutrinos
0209   std::vector<fastjet::PseudoJet> fjJetInputs;
0210   fjJetInputs.reserve(nStables);
0211   for (size_t i = 0, n = finalStateHandle->size(); i < n; ++i) {
0212     const reco::Candidate& p = finalStateHandle->at(i);
0213     if (p.status() != 1)
0214       continue;
0215     if (std::isnan(p.pt()) or p.pt() <= 0)
0216       continue;
0217 
0218     const int absId = std::abs(p.pdgId());
0219     if (absId == 12 or absId == 14 or absId == 16)
0220       continue;
0221     if (lepDauIdxs.find(i) != lepDauIdxs.end())
0222       continue;
0223 
0224     fjJetInputs.push_back(fastjet::PseudoJet(p.px(), p.py(), p.pz(), p.energy()));
0225     fjJetInputs.back().set_user_index(i);
0226   }
0227   //// Also don't forget to put B hadrons
0228   for (auto index : bHadronIdxs) {
0229     const reco::Candidate& p = genParticleHandle->at(abs(index + 1));
0230     if (std::isnan(p.pt()) or p.pt() <= 0)
0231       continue;
0232 
0233     const double scale = 1e-20 / p.p();
0234     fjJetInputs.push_back(fastjet::PseudoJet(p.px() * scale, p.py() * scale, p.pz() * scale, p.energy() * scale));
0235     fjJetInputs.back().set_user_index(index);
0236   }
0237 
0238   //// Run the jet algorithm
0239   fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *fjJetDef_);
0240   std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(minJetPt_));
0241 
0242   /// Build jets
0243   jets->reserve(fjJets.size());
0244   std::vector<size_t> bjetIdxs, ljetIdxs;
0245   for (auto& fjJet : fjJets) {
0246     if (abs(fjJet.eta()) > maxJetEta_)
0247       continue;
0248 
0249     // Get jet constituents from fastJet
0250     const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
0251     // Convert to CandidatePtr
0252     std::vector<reco::CandidatePtr> constituents;
0253     bool hasBHadron = false;
0254     for (size_t j = 0, m = fjConstituents.size(); j < m; ++j) {
0255       const int index = fjConstituents[j].user_index();
0256       if (bHadronIdxs.find(index) != bHadronIdxs.end()) {
0257         hasBHadron = true;
0258         continue;
0259       }
0260       reco::CandidatePtr cand = finalStateHandle->ptrAt(index);
0261       constituents.push_back(cand);
0262     }
0263 
0264     const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
0265     reco::GenJet genJet;
0266     reco::writeSpecific(genJet, jetP4, genVertex_, constituents);
0267 
0268     const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
0269     genJet.setJetArea(jetArea);
0270     if (hasBHadron) {
0271       genJet.setPdgId(5);
0272       bjetIdxs.push_back(jets->size());
0273     } else {
0274       ljetIdxs.push_back(jets->size());
0275     }
0276 
0277     jets->push_back(genJet);
0278   }
0279 
0280   // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
0281   // NOTE : A C++ trick, use do-while instead of long-nested if-statements.
0282   do {
0283     if (bjetIdxs.size() < 2)
0284       break;
0285 
0286     // Note : we will do dilepton or semilepton channel only
0287     if (leptons->size() == 2 and neutrinos->size() >= 2) {
0288       // Start from dilepton channel
0289       const int q1 = leptons->at(0).charge();
0290       const int q2 = leptons->at(1).charge();
0291       if (q1 * q2 > 0)
0292         break;
0293 
0294       const auto& lepton1 = q1 > 0 ? leptons->at(0) : leptons->at(1);
0295       const auto& lepton2 = q1 > 0 ? leptons->at(1) : leptons->at(0);
0296 
0297       if (lepton1.pt() < minLeptonPtDilepton_ or std::abs(lepton1.eta()) > maxLeptonEtaDilepton_)
0298         break;
0299       if (lepton2.pt() < minLeptonPtDilepton_ or std::abs(lepton2.eta()) > maxLeptonEtaDilepton_)
0300         break;
0301       if ((lepton1.p4() + lepton2.p4()).mass() < minDileptonMassDilepton_)
0302         break;
0303 
0304       double dm = 1e9;
0305       int selNu1 = -1, selNu2 = -1;
0306       for (int i = 0; i < 2; ++i) {  // Consider leading 2 neutrinos. Neutrino vector is already sorted by pT
0307         const double dm1 = std::abs((lepton1.p4() + neutrinos->at(i).p4()).mass() - wMass_);
0308         for (int j = 0; j < 2; ++j) {
0309           if (i == j)
0310             continue;
0311           const double dm2 = std::abs((lepton2.p4() + neutrinos->at(j).p4()).mass() - wMass_);
0312           const double newDm = dm1 + dm2;
0313 
0314           if (newDm < dm) {
0315             dm = newDm;
0316             selNu1 = i;
0317             selNu2 = j;
0318           }
0319         }
0320       }
0321       if (dm >= 1e9)
0322         break;
0323 
0324       const auto& nu1 = neutrinos->at(selNu1);
0325       const auto& nu2 = neutrinos->at(selNu2);
0326       const auto w1LVec = lepton1.p4() + nu1.p4();
0327       const auto w2LVec = lepton2.p4() + nu2.p4();
0328 
0329       // Contiue to top quarks
0330       dm = 1e9;  // Reset once again for top combination.
0331       int selB1 = -1, selB2 = -1;
0332       for (auto i : bjetIdxs) {
0333         const double dm1 = std::abs((w1LVec + jets->at(i).p4()).mass() - tMass_);
0334         for (auto j : bjetIdxs) {
0335           if (i == j)
0336             continue;
0337           const double dm2 = std::abs((w2LVec + jets->at(j).p4()).mass() - tMass_);
0338           const double newDm = dm1 + dm2;
0339 
0340           if (newDm < dm) {
0341             dm = newDm;
0342             selB1 = i;
0343             selB2 = j;
0344           }
0345         }
0346       }
0347       if (dm >= 1e9)
0348         break;
0349 
0350       const auto& bJet1 = jets->at(selB1);
0351       const auto& bJet2 = jets->at(selB2);
0352       const auto t1LVec = w1LVec + bJet1.p4();
0353       const auto t2LVec = w2LVec + bJet2.p4();
0354 
0355       // Recalculate lepton charges (needed due to initialization of leptons wrt sign of a charge)
0356       const int lepQ1 = lepton1.charge();
0357       const int lepQ2 = lepton2.charge();
0358 
0359       // Put all of them into candidate collection
0360       reco::GenParticle t1(lepQ1 * 2 / 3., t1LVec, genVertex_, lepQ1 * 6, 3, false);
0361       reco::GenParticle w1(lepQ1, w1LVec, genVertex_, lepQ1 * 24, 3, true);
0362       reco::GenParticle b1(-lepQ1 / 3., bJet1.p4(), genVertex_, lepQ1 * 5, 1, true);
0363       reco::GenParticle l1(lepQ1, lepton1.p4(), genVertex_, lepton1.pdgId(), 1, true);
0364       reco::GenParticle n1(0, nu1.p4(), genVertex_, nu1.pdgId(), 1, true);
0365 
0366       reco::GenParticle t2(lepQ2 * 2 / 3., t2LVec, genVertex_, lepQ2 * 6, 3, false);
0367       reco::GenParticle w2(lepQ2, w2LVec, genVertex_, lepQ2 * 24, 3, true);
0368       reco::GenParticle b2(-lepQ2 / 3., bJet2.p4(), genVertex_, lepQ2 * 5, 1, true);
0369       reco::GenParticle l2(lepQ2, lepton2.p4(), genVertex_, lepton2.pdgId(), 1, true);
0370       reco::GenParticle n2(0, nu2.p4(), genVertex_, nu2.pdgId(), 1, true);
0371 
0372       pseudoTop->push_back(t1);
0373       pseudoTop->push_back(t2);
0374 
0375       pseudoTop->push_back(w1);
0376       pseudoTop->push_back(b1);
0377       pseudoTop->push_back(l1);
0378       pseudoTop->push_back(n1);
0379 
0380       pseudoTop->push_back(w2);
0381       pseudoTop->push_back(b2);
0382       pseudoTop->push_back(l2);
0383       pseudoTop->push_back(n2);
0384     } else if (!leptons->empty() and !neutrinos->empty() and ljetIdxs.size() >= 2) {
0385       // Then continue to the semi-leptonic channel
0386       const auto& lepton = leptons->at(0);
0387       if (lepton.pt() < minLeptonPtSemilepton_ or std::abs(lepton.eta()) > maxLeptonEtaSemilepton_)
0388         break;
0389 
0390       // Skip event if there are veto leptons
0391       bool hasVetoLepton = false;
0392       for (auto vetoLeptonItr = next(leptons->begin()); vetoLeptonItr != leptons->end(); ++vetoLeptonItr) {
0393         if (vetoLeptonItr->pt() > minVetoLeptonPtSemilepton_ and
0394             std::abs(vetoLeptonItr->eta()) < maxVetoLeptonEtaSemilepton_) {
0395           hasVetoLepton = true;
0396         }
0397       }
0398       if (hasVetoLepton)
0399         break;
0400 
0401       // Calculate MET
0402       double metX = 0, metY = 0;
0403       for (const auto& neutrino : *neutrinos) {
0404         metX += neutrino.px();
0405         metY += neutrino.py();
0406       }
0407       const double metPt = std::hypot(metX, metY);
0408       if (metPt < minMETSemiLepton_)
0409         break;
0410 
0411       const double mtW = std::sqrt(2 * lepton.pt() * metPt * cos(lepton.phi() - atan2(metX, metY)));
0412       if (mtW < minMtWSemiLepton_)
0413         break;
0414 
0415       // Solve pz to give wMass_^2 = (lepE+energy)^2 - (lepPx+metX)^2 - (lepPy+metY)^2 - (lepPz+pz)^2
0416       // -> (wMass_^2)/2 = lepE*sqrt(metPt^2+pz^2) - lepPx*metX - lepPy*metY - lepPz*pz
0417       // -> C := (wMass_^2)/2 + lepPx*metX + lepPy*metY
0418       // -> lepE^2*(metPt^2+pz^2) = C^2 + 2*C*lepPz*pz + lepPz^2*pz^2
0419       // -> (lepE^2-lepPz^2)*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
0420       // -> lepPt^2*pz^2 - 2*C*lepPz*pz - C^2 + lepE^2*metPt^2 = 0
0421       const double lpz = lepton.pz(), le = lepton.energy(), lpt = lepton.pt();
0422       const double cf = (wMass_ * wMass_) / 2 + lepton.px() * metX + lepton.py() * metY;
0423       const double det = cf * cf * lpz * lpz - lpt * lpt * (le * le * metPt * metPt - cf * cf);
0424       const double pz = (cf * lpz + (cf < 0 ? -sqrt(det) : sqrt(det))) / lpt / lpt;
0425       const reco::Candidate::LorentzVector nu1P4(metX, metY, pz, std::hypot(metPt, pz));
0426       const auto w1LVec = lepton.p4() + nu1P4;
0427 
0428       // Continue to build leptonic pseudo top
0429       double minDR = 1e9;
0430       int selB1 = -1;
0431       for (auto b1Itr = bjetIdxs.begin(); b1Itr != bjetIdxs.end(); ++b1Itr) {
0432         const double dR = deltaR(jets->at(*b1Itr).p4(), w1LVec);
0433         if (dR < minDR) {
0434           selB1 = *b1Itr;
0435           minDR = dR;
0436         }
0437       }
0438       if (selB1 == -1)
0439         break;
0440       const auto& bJet1 = jets->at(selB1);
0441       const auto t1LVec = w1LVec + bJet1.p4();
0442 
0443       // Build hadronic pseudo W, take leading two jets (ljetIdxs are (implicitly) sorted by pT)
0444       const auto& wJet1 = jets->at(ljetIdxs[0]);
0445       const auto& wJet2 = jets->at(ljetIdxs[1]);
0446       const auto w2LVec = wJet1.p4() + wJet2.p4();
0447 
0448       // Contiue to top quarks
0449       double dm = 1e9;
0450       int selB2 = -1;
0451       for (auto i : bjetIdxs) {
0452         if (int(i) == selB1)
0453           continue;
0454         const double newDm = std::abs((w2LVec + jets->at(i).p4()).mass() - tMass_);
0455         if (newDm < dm) {
0456           dm = newDm;
0457           selB2 = i;
0458         }
0459       }
0460       if (dm >= 1e9)
0461         break;
0462 
0463       const auto& bJet2 = jets->at(selB2);
0464       const auto t2LVec = w2LVec + bJet2.p4();
0465 
0466       const int q = lepton.charge();
0467       // Put all of them into candidate collection
0468       reco::GenParticle t1(q * 2 / 3., t1LVec, genVertex_, q * 6, 3, false);
0469       reco::GenParticle w1(q, w1LVec, genVertex_, q * 24, 3, true);
0470       reco::GenParticle b1(-q / 3., bJet1.p4(), genVertex_, q * 5, 1, true);
0471       reco::GenParticle l1(q, lepton.p4(), genVertex_, lepton.pdgId(), 1, true);
0472       reco::GenParticle n1(0, nu1P4, genVertex_, q * (std::abs(lepton.pdgId()) + 1), 1, true);
0473 
0474       reco::GenParticle t2(-q * 2 / 3., t2LVec, genVertex_, -q * 6, 3, false);
0475       reco::GenParticle w2(-q, w2LVec, genVertex_, -q * 24, 3, true);
0476       reco::GenParticle b2(0, bJet2.p4(), genVertex_, -q * 5, 1, true);
0477       reco::GenParticle u2(0, wJet1.p4(), genVertex_, -2 * q, 1, true);
0478       reco::GenParticle d2(0, wJet2.p4(), genVertex_, q, 1, true);
0479 
0480       pseudoTop->push_back(t1);
0481       pseudoTop->push_back(t2);
0482 
0483       pseudoTop->push_back(w1);
0484       pseudoTop->push_back(b1);
0485       pseudoTop->push_back(l1);
0486       pseudoTop->push_back(n1);
0487 
0488       pseudoTop->push_back(w2);
0489       pseudoTop->push_back(b2);
0490       pseudoTop->push_back(u2);
0491       pseudoTop->push_back(d2);
0492     }
0493   } while (false);
0494 
0495   if (pseudoTop->size() == 10)  // If pseudtop decay tree is completed
0496   {
0497     // t->W+b
0498     pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 2));  // t->W
0499     pseudoTop->at(0).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 3));  // t->b
0500     pseudoTop->at(2).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0));    // t->W
0501     pseudoTop->at(3).addMother(reco::GenParticleRef(pseudoTopRefHandle, 0));    // t->b
0502 
0503     // W->lv or W->jj
0504     pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 4));
0505     pseudoTop->at(2).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 5));
0506     pseudoTop->at(4).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
0507     pseudoTop->at(5).addMother(reco::GenParticleRef(pseudoTopRefHandle, 2));
0508 
0509     // tbar->W-b
0510     pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 6));
0511     pseudoTop->at(1).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 7));
0512     pseudoTop->at(6).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
0513     pseudoTop->at(7).addMother(reco::GenParticleRef(pseudoTopRefHandle, 1));
0514 
0515     // W->jj
0516     pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 8));
0517     pseudoTop->at(6).addDaughter(reco::GenParticleRef(pseudoTopRefHandle, 9));
0518     pseudoTop->at(8).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
0519     pseudoTop->at(9).addMother(reco::GenParticleRef(pseudoTopRefHandle, 6));
0520   }
0521 
0522   event.put(std::move(neutrinos), "neutrinos");
0523   event.put(std::move(leptons), "leptons");
0524   event.put(std::move(jets), "jets");
0525 
0526   event.put(std::move(pseudoTop));
0527 }
0528 
0529 bool PseudoTopProducer::isFromHadron(const reco::Candidate* p) const {
0530   for (size_t i = 0, n = p->numberOfMothers(); i < n; ++i) {
0531     const reco::Candidate* mother = p->mother(i);
0532     if (!mother)
0533       continue;
0534     if (mother->status() == 4 or mother->numberOfMothers() == 0)
0535       continue;  // Skip incident beam
0536     const int pdgId = abs(mother->pdgId());
0537 
0538     if (pdgId > 100)
0539       return true;
0540     else if (isFromHadron(mother))
0541       return true;
0542   }
0543   return false;
0544 }
0545 
0546 bool PseudoTopProducer::isBHadron(const reco::Candidate* p) const {
0547   const unsigned int absPdgId = abs(p->pdgId());
0548   if (!isBHadron(absPdgId))
0549     return false;
0550 
0551   // Do not consider this particle if it has B hadron daughter
0552   // For example, B* -> B0 + photon; then we drop B* and take B0 only
0553   for (int i = 0, n = p->numberOfDaughters(); i < n; ++i) {
0554     const reco::Candidate* dau = p->daughter(i);
0555     if (isBHadron(abs(dau->pdgId())))
0556       return false;
0557   }
0558 
0559   return true;
0560 }
0561 
0562 bool PseudoTopProducer::isBHadron(const unsigned int absPdgId) const {
0563   if (absPdgId <= 100)
0564     return false;  // Fundamental particles and MC internals
0565   if (absPdgId >= 1000000000)
0566     return false;  // Nuclei, +-10LZZZAAAI
0567 
0568   // General form of PDG ID is 7 digit form
0569   // +- n nr nL nq1 nq2 nq3 nJ
0570   //const int nJ = absPdgId % 10; // Spin
0571   const int nq3 = (absPdgId / 10) % 10;
0572   const int nq2 = (absPdgId / 100) % 10;
0573   const int nq1 = (absPdgId / 1000) % 10;
0574 
0575   if (nq3 == 0)
0576     return false;  // Diquarks
0577   if (nq1 == 0 and nq2 == 5)
0578     return true;  // B mesons
0579   if (nq1 == 5)
0580     return true;  // B baryons
0581 
0582   return false;
0583 }
0584 
0585 #include "FWCore/Framework/interface/MakerMacros.h"
0586 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!