// -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/PartonicTops.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// Helper projection for struct PseudoBoostedTop : public FinalState { // Default constructor PseudoBoostedTop (double lepMinPt = 45., double lepMaxEta = 2.1, double jetMinPt = 30., double jetMaxEta = 2.4, double bjetR = 0.5, double tjetR = 0.8) : FinalState(), _lepMinPt(lepMinPt), _lepMaxEta(lepMaxEta), _jetMinPt(jetMinPt), _jetMaxEta(jetMaxEta), _bjetR(bjetR), _tjetR(tjetR) { setName("PseudoBoostedTop"); } enum DecayMode { CH_HADRON, CH_MUON, CH_ELECTRON, CH_TAU }; enum TTbarMode { CH_DILEP, CH_ALLHADRON, CH_SEMILEP_MU, CH_SEMILEP_EL, CH_SEMILEP_TAU }; TTbarMode mode() const { if ((_topDecay == CH_MUON && _antitopDecay == CH_HADRON) || (_topDecay == CH_HADRON && _antitopDecay == CH_MUON)) { return CH_SEMILEP_MU; } else if ((_topDecay == CH_ELECTRON && _antitopDecay == CH_HADRON) || (_topDecay == CH_HADRON && _antitopDecay == CH_ELECTRON)) { return CH_SEMILEP_EL; } else if ((_topDecay == CH_TAU && _antitopDecay == CH_HADRON) || (_topDecay == CH_HADRON && _antitopDecay == CH_TAU)) { return CH_SEMILEP_TAU; } else if (_topDecay == CH_HADRON && _antitopDecay == CH_HADRON) { return CH_ALLHADRON; } else { return CH_DILEP; } } Particle partonTop() const { if (_antitopDecay == CH_HADRON) return _antitop; else if (_topDecay == CH_HADRON) return _top; else return nullptr; } Jet particleTop() const { return _particleTop; } bool isGenTopJet() const { return _isGenTopJet; } bool passParticle() const { return _passParticle; } // Determine if lepton comes from hard top decay bool isFromTop(const GenParticle* p) { bool fromTop = false; const int pdgId = p->pdg_id(); GenVertex* prodVtx = p->production_vertex(); if (prodVtx != nullptr) { for (const GenParticle* p1 : Rivet::particles(prodVtx, HepMC::parents)) { const int pdgIdParent = p1->pdg_id(); if (pdgIdParent == 24) { // parent is a W+ --> leptonic top if (pdgId == -11) _topDecay = CH_ELECTRON; if (pdgId == -13) _topDecay = CH_MUON; if (pdgId == -15) _topDecay = CH_TAU; fromTop = true; } else if (pdgIdParent == -24) { // parent is a W- --> leptonic antitop if (pdgId == 11) _antitopDecay = CH_ELECTRON; if (pdgId == 13) _antitopDecay = CH_MUON; if (pdgId == 15) _antitopDecay = CH_TAU; fromTop = true; } } } return fromTop; } bool fromHardMuon(const GenParticle* p) { GenVertex* prodVtx = p->production_vertex(); if (prodVtx == nullptr) return false; for (const GenParticle* ancestor : Rivet::particles(prodVtx, HepMC::ancestors)) { if (std::abs(ancestor->pdg_id()) == PID::MUON && ancestor->status() == 3) ///< @todo Unsafe!!! return true; } return false; } const double _lepMinPt, _lepMaxEta; const double _jetMinPt, _jetMaxEta; const double _bjetR, _tjetR; bool _isValid; DecayMode _topDecay, _antitopDecay; Particle _top, _antitop; Jet _particleTop; bool _isParticleLep, _isGenBJet, _isGenTopJet, _passParticle; double _mtt; }; /// @brief Boosted ttbar in pp collisions at sqrt{s} = 8 TeV /// /// @warning Relies on parton quantities and unsafe MC status codes, and weights are not used properly /// /// In this analysis, the parton-level phase space is defined using the following (UNSAFE!) object definitions: /// * Lepton: An electron, muon, or tau (PDG ID = $\pm11$, $\pm13$, or $\pm15$) originating from the /// decay of a W boson (PDG ID of parent particle is $\pm24$). /// * Top quark: A top quark (PDG ID = $\pm6$) with status 3 or 22. /// /// The parton-level phase space is defined by requiring exactly one lepton, which must be an electron or muon. /// For ${\rm t\bar{t}}$ Monte Carlo, this is equivalent to requiring the event to be semileptonic at parton level. /// The lepton charge is used to determine whether it comes from the decay of a top or antitop. The top of /// opposite sign (which by construct must have decayed hadronically) is defined as the parton-level top quark. /// The parton-level top quark is required to have $p_{T} > 400$~GeV. /// /// In this analysis, the particle-level phase space is defined using the following object definitions: /// * Lepton: An electron or muon (PDG ID = $\pm11$ or $\pm13$) originating from the decay of a W boson /// (PDG ID of parent particle is $\pm24$). The lepton is required to have $p_{T} > 45$~GeV and $|\eta| < 2.1$. /// * B Jet Candidate: Gen AK5 jets are formed by clustering the final state particles in the event using /// the anti-$k_{T}$ algorithm with distance parameter $R=0.5$. Here, 'final state particles' refers to /// all status 1 particles in the event, excepting neutrinos and particles originating from the decay /// of the muon defined above. The gen AK5 jet is required to have $p_{T} > 30$~GeV and $|\eta| < 2.4$. /// Gen AK5 jets in the same hemisphere as the lepton ($\Delta {\rm R(e/\mu,~jet)} < \pi/2$) are /// defined as b-jet candidates. /// * Top Jet Candidate: Gen CA8 jets are formed by clustering the final state particles in the event using /// the Cambridge-Aachen algorithm with distance parameter $R=0.8$. Here, 'final state particles' refers /// to all status 1 particles in the event, excepting neutrinos and particles originating from the decay /// of the muon defined above. The gen CA8 jet is required to have $p_{T} > 30$~GeV and $|\eta| < 2.4$. /// Gen CA8 jets which have $p_{T} > 400$~GeV, 140 GeV $<$ mass $<$ 250 GeV, and are in the opposite /// hemisphere from the lepton ($\Delta{\rm R(e/\mu,~jet)} > \pi/2$) are defined as top jet candidates. /// /// The particle-level phase space is defined by requiring $\geq 1$ b jet candidate, $\geq 1$ top jet candidate, /// and exactly one lepton. This is in addition to the parton-level semileptonic requirement. The highest-$p_{T}$ /// top jet candidate is defined as the particle-level t jet. /// class CMS_2015_I1388555 : public Analysis { public: // Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2015_I1388555); // Book histograms and initialise projections before the run void init() { declare(PseudoBoostedTop(45., 2.1, 30., 2.4, 0.5, 0.8), "ttbar"); _hEl_topPt_parton = bookHisto1D("d01-x01-y01"); // dsigma/dpt(top) _hEl_topPt_particle = bookHisto1D("d02-x01-y01"); // dsigma/dpt(top) _hEl_topY_parton = bookHisto1D("d03-x01-y01"); // dsigma/dy(top) _hEl_topY_particle = bookHisto1D("d04-x01-y01"); // dsigma/dy(top) _hMu_topPt_parton = bookHisto1D("d05-x01-y01"); // dsigma/dpt(top) _hMu_topPt_particle = bookHisto1D("d06-x01-y01"); // dsigma/dpt(top) _hMu_topY_parton = bookHisto1D("d07-x01-y01"); // dsigma/dy(top) _hMu_topY_particle = bookHisto1D("d08-x01-y01"); // dsigma/dy(top) _hComb_topPt_parton = bookHisto1D("d09-x01-y01"); // dsigma/dpt(top) _hComb_topPt_particle = bookHisto1D("d10-x01-y01"); // dsigma/dpt(top) _hComb_topY_parton = bookHisto1D("d11-x01-y01"); // dsigma/dy(top) _hComb_topY_particle = bookHisto1D("d12-x01-y01"); // dsigma/dy(top) _hEl_topPt_parton_norm = bookHisto1D("d13-x01-y01"); // 1/sigma dsigma/dpt(top) _hEl_topPt_particle_norm = bookHisto1D("d14-x01-y01"); // 1/sigma dsigma/dpt(top) _hEl_topY_parton_norm = bookHisto1D("d15-x01-y01"); // 1/sigma dsigma/dy(top) _hEl_topY_particle_norm = bookHisto1D("d16-x01-y01"); // 1/sigma dsigma/dy(top) _hMu_topPt_parton_norm = bookHisto1D("d17-x01-y01"); // 1/sigma dsigma/dpt(top) _hMu_topPt_particle_norm = bookHisto1D("d18-x01-y01"); // 1/sigma dsigma/dpt(top) _hMu_topY_parton_norm = bookHisto1D("d19-x01-y01"); // 1/sigma dsigma/dy(top) _hMu_topY_particle_norm = bookHisto1D("d20-x01-y01"); // 1/sigma dsigma/dy(top) _hComb_topPt_parton_norm = bookHisto1D("d21-x01-y01"); // 1/sigma dsigma/dpt(top) _hComb_topPt_particle_norm = bookHisto1D("d22-x01-y01"); // 1/sigma dsigma/dpt(top) _hComb_topY_parton_norm = bookHisto1D("d23-x01-y01"); // 1/sigma dsigma/dy(top) _hComb_topY_particle_norm = bookHisto1D("d24-x01-y01"); // 1/sigma dsigma/dy(top) nMu = 0, nEl = 0; nPassParton_mu = 0, nPassParticle_mu = 0; nPassParton_el = 0, nPassParticle_el = 0; } /// Perform the per-event analysis void analyze(const Event& event) { // Get the ttbar candidate const PseudoBoostedTop& ttbar = applyProjection(event, "ttbar"); if (ttbar.mode() == PseudoBoostedTop::CH_DILEP || ttbar.mode() == PseudoBoostedTop::CH_ALLHADRON || ttbar.mode() == PseudoBoostedTop::CH_SEMILEP_TAU) vetoEvent; const double weight = event.weight(); // Do selection const FourMomentum& partonTopP4 = ttbar.partonTop().momentum(); if (ttbar.mode() == PseudoBoostedTop::CH_SEMILEP_MU) { // pass loose parton, muon channel nMu += 1; _hMu_topPt_parton->fill(partonTopP4.pT()/GeV, weight); _hMu_topPt_parton_norm->fill(partonTopP4.pT()/GeV, weight); _hComb_topPt_parton->fill(partonTopP4.pT()/GeV, weight); _hComb_topPt_parton_norm->fill(partonTopP4.pT()/GeV, weight); if (partonTopP4.pT() > 400*GeV) { nPassParton_mu += 1; _hMu_topY_parton->fill(partonTopP4.rapidity(), weight); _hMu_topY_parton_norm->fill(partonTopP4.rapidity(), weight); _hComb_topY_parton->fill(partonTopP4.rapidity(), weight); _hComb_topY_parton_norm->fill(partonTopP4.rapidity(), weight); } if (ttbar.isGenTopJet()) { const FourMomentum& particleTopP4 = ttbar.particleTop().momentum(); _hMu_topPt_particle->fill(particleTopP4.pT()/GeV, weight); _hMu_topPt_particle_norm->fill(particleTopP4.pT()/GeV, weight); _hComb_topPt_particle->fill(particleTopP4.pT()/GeV, weight); _hComb_topPt_particle_norm->fill(particleTopP4.pT()/GeV, weight); if (ttbar.passParticle()){ nPassParticle_mu += 1; _hMu_topY_particle->fill(particleTopP4.rapidity(), weight); _hMu_topY_particle_norm->fill(particleTopP4.rapidity(), weight); _hComb_topY_particle->fill(particleTopP4.rapidity(), weight); _hComb_topY_particle_norm->fill(particleTopP4.rapidity(), weight); } } } else if (ttbar.mode() == PseudoBoostedTop::CH_SEMILEP_EL) { nEl += 1; _hEl_topPt_parton->fill(partonTopP4.pT()/GeV, weight); _hEl_topPt_parton_norm->fill(partonTopP4.pT()/GeV, weight); _hComb_topPt_parton->fill(partonTopP4.pT()/GeV, weight); _hComb_topPt_parton_norm->fill(partonTopP4.pT()/GeV, weight); if (partonTopP4.pT() >= 400*GeV) { nPassParton_el += 1; _hEl_topY_parton->fill(partonTopP4.rapidity(), weight); _hEl_topY_parton_norm->fill(partonTopP4.rapidity(), weight); _hComb_topY_parton->fill(partonTopP4.rapidity(), weight); _hComb_topY_parton_norm->fill(partonTopP4.rapidity(), weight); } if (ttbar.isGenTopJet()) { const FourMomentum& particleTopP4 = ttbar.particleTop().momentum(); _hEl_topPt_particle->fill(particleTopP4.pT()/GeV, weight); _hEl_topPt_particle_norm->fill(particleTopP4.pT()/GeV, weight); _hComb_topPt_particle->fill(particleTopP4.pT()/GeV, weight); _hComb_topPt_particle_norm->fill(particleTopP4.pT()/GeV, weight); if (ttbar.passParticle()) { nPassParticle_el += 1; _hEl_topY_particle->fill(particleTopP4.rapidity(), weight); _hEl_topY_particle_norm->fill(particleTopP4.rapidity(), weight); _hComb_topY_particle->fill(particleTopP4.rapidity(), weight); _hComb_topY_particle_norm->fill(particleTopP4.rapidity(), weight); } } } } // Normalise histograms etc., after the run void finalize() { /// @todo What are the normalisation factors? What about generator xsec? const double xs_mu_parton = 252.89 * 1000. * nPassParton_mu / nMu; const double xs_mu_particle = 252.89 * 1000. * nPassParticle_mu / nMu; const double xs_el_parton = 252.89 * 1000. * nPassParton_el / nEl; const double xs_el_particle = 252.89 * 1000. * nPassParticle_el / nEl; const double xs_comb_parton = 252.89 * 1000. * (nPassParton_el + nPassParton_mu) / (nEl + nMu); const double xs_comb_particle = 252.89 * 1000. * (nPassParticle_el + nPassParton_mu) / (nEl + nMu); normalize(_hMu_topPt_parton, xs_mu_parton, false); normalize(_hMu_topPt_particle, xs_mu_particle, false); normalize(_hMu_topY_parton, xs_mu_parton, false); normalize(_hMu_topY_particle, xs_mu_particle, false); normalize(_hEl_topPt_parton, xs_el_parton, false); normalize(_hEl_topPt_particle, xs_el_particle, false); normalize(_hEl_topY_parton, xs_el_parton, false); normalize(_hEl_topY_particle, xs_el_particle, false); normalize(_hComb_topPt_parton, xs_comb_parton, false); normalize(_hComb_topPt_particle, xs_comb_particle, false); normalize(_hComb_topY_parton, xs_comb_parton, false); normalize(_hComb_topY_particle, xs_comb_particle, false); normalize(_hMu_topPt_parton_norm, 1.0, false); normalize(_hMu_topPt_particle_norm, 1.0, false); normalize(_hMu_topY_parton_norm, 1.0, false); normalize(_hMu_topY_particle_norm, 1.0, false); normalize(_hEl_topPt_parton_norm, 1.0, false); normalize(_hEl_topPt_particle_norm, 1.0, false); normalize(_hEl_topY_parton_norm, 1.0, false); normalize(_hEl_topY_particle_norm, 1.0, false); normalize(_hComb_topPt_parton_norm, 1.0, false); normalize(_hComb_topPt_particle_norm, 1.0, false); normalize(_hComb_topY_parton_norm, 1.0, false); normalize(_hComb_topY_particle_norm, 1.0, false); } private: Histo1DPtr _hMu_topPt_parton; Histo1DPtr _hMu_topPt_particle; Histo1DPtr _hMu_topY_parton; Histo1DPtr _hMu_topY_particle; Histo1DPtr _hEl_topPt_parton; Histo1DPtr _hEl_topPt_particle; Histo1DPtr _hEl_topY_parton; Histo1DPtr _hEl_topY_particle; Histo1DPtr _hComb_topPt_parton; Histo1DPtr _hComb_topPt_particle; Histo1DPtr _hComb_topY_parton; Histo1DPtr _hComb_topY_particle; Histo1DPtr _hMu_topPt_parton_norm; Histo1DPtr _hMu_topPt_particle_norm; Histo1DPtr _hMu_topY_parton_norm; Histo1DPtr _hMu_topY_particle_norm; Histo1DPtr _hEl_topPt_parton_norm; Histo1DPtr _hEl_topPt_particle_norm; Histo1DPtr _hEl_topY_parton_norm; Histo1DPtr _hEl_topY_particle_norm; Histo1DPtr _hComb_topPt_parton_norm; Histo1DPtr _hComb_topPt_particle_norm; Histo1DPtr _hComb_topY_parton_norm; Histo1DPtr _hComb_topY_particle_norm; /// @todo These need to be weights int nMu, nEl, nPassParton_mu, nPassParticle_mu, nPassParton_el, nPassParticle_el; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2015_I1388555); }