Rivet analyses
Differential top quark pair production cross-sections in pp collisions at $\sqrt{s} = 8$ TeV
Experiment: CMS (LHC)
Inspire ID: 1370682
Status: VALIDATED
Authors: - Javier Fernandez - Jungwan John Goh - Efe Yazgan - Markus Seidel - James Keaveney - Elvire Bouvier - Benedikt Maier
References: - Eur.Phys.J. C75 (2015) 542 (the analysis) - Expt page: CMS-PAS-TOP-15-011 (Rivet implementation)
Beams: p+ p+
Beam energies: (4000.0, 4000.0)GeV
Run details: - pp QCD interactions at $\sqrt{s} = 8$ TeV. Data collected by CMS during the year 2012.
The reconstruction of particle-level top quarks and anti-quarks is implemented in this module. Measurements at $\sqrt{s} = 8~\TeV$ are based on parton-level information in the full phase space using MADGRAPH+PYTHIA6. To match the particle-level top quark distributions to the measurements unfolded to the parton-level, a correction function to the particle-level distributions, derived using the same MADGRAPH+PYTHIA6 configuration that was used for the original measurement of the data points, is applied. Using the same MC configuration as used for the unfolding to correct back the parton-level to particle-level, the model dependence introduced in unfolding to parton-level and extrapolating the measurement to the full phase space is eliminated. See the paper for full object selection and correction details.
Source
code:CMS_2015_I1370682.cc
#include "Rivet/Analysis.hh"
#include "Rivet/Math/LorentzTrans.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
namespace Rivet {
namespace { //< only visible in this compilation unit
/// @brief Pseudo top finder
///
/// Find top quark in the particle level.
/// The definition is based on the agreement at the LHC working group.
class PseudoTop : public FinalState {
public:
/// @name Standard constructors and destructors.
/// @{
/// The default constructor. May specify the minimum and maximum
/// pseudorapidity \f$ \eta \f$ and the min \f$ p_T \f$ (in GeV).
PseudoTop(double lepR = 0.1,
double lepMinPt = 20,
double lepMaxEta = 2.4,
double jetR = 0.4,
double jetMinPt = 30,
double jetMaxEta = 4.7)
: FinalState(),
_lepR(lepR),
_lepMinPt(lepMinPt),
_lepMaxEta(lepMaxEta),
_jetR(jetR),
_jetMinPt(jetMinPt),
_jetMaxEta(jetMaxEta) {
setName("PseudoTop");
}
enum TTbarMode { CH_NONE = -1, CH_FULLHADRON = 0, CH_SEMILEPTON, CH_FULLLEPTON };
enum DecayMode { CH_HADRON = 0, CH_MUON, CH_ELECTRON };
TTbarMode mode() const {
if (!_isValid) return CH_NONE;
if (_mode1 == CH_HADRON && _mode2 == CH_HADRON)
return CH_FULLHADRON;
else if (_mode1 != CH_HADRON && _mode2 != CH_HADRON)
return CH_FULLLEPTON;
else
return CH_SEMILEPTON;
}
virtual const DecayMode& mode1() const {
return _mode1;
}
virtual const DecayMode& mode2() const {
return _mode2;
}
/// Clone on the heap.
RIVET_DEFAULT_PROJ_CLONE(PseudoTop);
/// @}
/// Import to avoid warnings about overload-hiding
using Projection::operator=;
public:
virtual const Particle& t1() const {
return _t1;
}
virtual const Particle& t2() const {
return _t2;
}
virtual const Particle& b1() const {
return _b1;
}
virtual const Particle& b2() const {
return _b2;
}
virtual const Particles& wDecays1() const {
return _wDecays1;
}
virtual const Particles& wDecays2() const {
return _wDecays2;
}
virtual const Jets& jets() const {
return _jets;
}
virtual const Jets& bjets() const {
return _bjets;
}
virtual const Jets& ljets() const {
return _ljets;
}
// Apply the projection to the event
void project(const Event& e); // override; ///< @todo Re-enable when C++11 allowed
void cleanup(std::map<double, std::pair<size_t, size_t>>& v, const bool doCrossCleanup = false) const;
CmpState compare(const Projection& p) const;
private:
const double _lepR, _lepMinPt, _lepMaxEta;
const double _jetR, _jetMinPt, _jetMaxEta;
//constexpr ///< @todo Re-enable when C++11 allowed
static double _tMass; // = 172.5*GeV; ///< @todo Re-enable when C++11 allowed
//constexpr ///< @todo Re-enable when C++11 allowed
static double _wMass; // = 80.4*GeV; ///< @todo Re-enable when C++11 allowed
private:
bool _isValid;
DecayMode _mode1, _mode2;
Particle _t1, _t2;
Particle _b1, _b2;
Particles _wDecays1, _wDecays2;
Jets _jets, _bjets, _ljets;
};
// More implementation below the analysis code
}
/// Pseudo-top analysis from CMS
class CMS_2015_I1370682 : public Analysis {
public:
CMS_2015_I1370682()
: Analysis("CMS_2015_I1370682"), _applyCorrection(true), _doShapeOnly(true) { }
void init() {
declare(PseudoTop(0.1, 20, 2.4, 0.5, 30, 2.4), "ttbar");
// Lepton + Jet channel
book(_hSL_topPt, "d15-x01-y01"); // 1/sigma dsigma/dpt(top)
book(_hSL_topPtTtbarSys, "d16-x01-y01"); // 1/sigma dsigma/dpt*(top)
book(_hSL_topY, "d17-x01-y01"); // 1/sigma dsigma/dy(top)
book(_hSL_ttbarDelPhi, "d18-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
book(_hSL_topPtLead, "d19-x01-y01"); // 1/sigma dsigma/dpt(t1)
book(_hSL_topPtSubLead, "d20-x01-y01"); // 1/sigma dsigma/dpt(t2)
book(_hSL_ttbarPt, "d21-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
book(_hSL_ttbarY, "d22-x01-y01"); // 1/sigma dsigma/dy(ttbar)
book(_hSL_ttbarMass, "d23-x01-y01"); // 1/sigma dsigma/dm(ttbar)
// Dilepton channel
book(_hDL_topPt, "d24-x01-y01"); // 1/sigma dsigma/dpt(top)
book(_hDL_topPtTtbarSys, "d25-x01-y01"); // 1/sigma dsigma/dpt*(top)
book(_hDL_topY, "d26-x01-y01"); // 1/sigma dsigma/dy(top)
book(_hDL_ttbarDelPhi, "d27-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
book(_hDL_topPtLead, "d28-x01-y01"); // 1/sigma dsigma/dpt(t1)
book(_hDL_topPtSubLead, "d29-x01-y01"); // 1/sigma dsigma/dpt(t2)
book(_hDL_ttbarPt, "d30-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
book(_hDL_ttbarY, "d31-x01-y01"); // 1/sigma dsigma/dy(ttbar)
book(_hDL_ttbarMass, "d32-x01-y01"); // 1/sigma dsigma/dm(ttbar)
}
void analyze(const Event& event) {
// Get the ttbar candidate
const PseudoTop& ttbar = apply<PseudoTop>(event, "ttbar");
if (ttbar.mode() == PseudoTop::CH_NONE) vetoEvent;
const FourMomentum& t1P4 = ttbar.t1().momentum();
const FourMomentum& t2P4 = ttbar.t2().momentum();
const double pt1 = std::max(t1P4.pT(), t2P4.pT());
const double pt2 = std::min(t1P4.pT(), t2P4.pT());
const double dPhi = deltaPhi(t1P4, t2P4);
const FourMomentum ttP4 = t1P4 + t2P4;
const FourMomentum t1P4AtCM = LorentzTransform::mkFrameTransformFromBeta(ttP4.betaVec())
.transform(t1P4);
if (ttbar.mode() == PseudoTop::CH_SEMILEPTON) {
const Particle lCand1 = ttbar.wDecays1()[0]; // w1 dau0 is the lepton in the PseudoTop
if (lCand1.pT() < 33 * GeV || lCand1.abseta() > 2.1) vetoEvent;
_hSL_topPt->fill(t1P4.pT());
_hSL_topPt->fill(t2P4.pT());
_hSL_topPtTtbarSys->fill(t1P4AtCM.pT());
_hSL_topY->fill(t1P4.rapidity());
_hSL_topY->fill(t2P4.rapidity());
_hSL_ttbarDelPhi->fill(dPhi);
_hSL_topPtLead->fill(pt1);
_hSL_topPtSubLead->fill(pt2);
_hSL_ttbarPt->fill(ttP4.pT());
_hSL_ttbarY->fill(ttP4.rapidity());
_hSL_ttbarMass->fill(ttP4.mass());
}
else if (ttbar.mode() == PseudoTop::CH_FULLLEPTON) {
const Particle lCand1 = ttbar.wDecays1()[0]; // dau0 are the lepton in the PseudoTop
const Particle lCand2 = ttbar.wDecays2()[0]; // dau0 are the lepton in the PseudoTop
if (lCand1.pT() < 20 * GeV || lCand1.abseta() > 2.4) vetoEvent;
if (lCand2.pT() < 20 * GeV || lCand2.abseta() > 2.4) vetoEvent;
_hDL_topPt->fill(t1P4.pT());
_hDL_topPt->fill(t2P4.pT());
_hDL_topPtTtbarSys->fill(t1P4AtCM.pT());
_hDL_topY->fill(t1P4.rapidity());
_hDL_topY->fill(t2P4.rapidity());
_hDL_ttbarDelPhi->fill(dPhi);
_hDL_topPtLead->fill(pt1);
_hDL_topPtSubLead->fill(pt2);
_hDL_ttbarPt->fill(ttP4.pT());
_hDL_ttbarY->fill(ttP4.rapidity());
_hDL_ttbarMass->fill(ttP4.mass());
}
}
void finalize() {
if (_applyCorrection) {
// Correction functions for TOP-12-028 paper, (parton bin height)/(pseudotop bin height)
const double ch15[] = {
5.473609, 4.941048, 4.173346, 3.391191, 2.785644, 2.371346, 2.194161, 2.197167,
};
const double ch16[] = {
5.470905, 4.948201, 4.081982, 3.225532, 2.617519, 2.239217, 2.127878, 2.185918,
};
const double ch17[] = {
10.003667, 4.546519, 3.828115, 3.601018, 3.522194,
3.524694, 3.600951, 3.808553, 4.531891, 9.995370,
};
const double ch18[] = {
4.406683,
4.054041,
3.885393,
4.213646,
};
const double ch19[] = {
6.182537, 5.257703, 4.422280, 3.568402, 2.889408, 2.415878, 2.189974, 2.173210,
};
const double ch20[] = {
5.199874, 4.693318, 3.902882, 3.143785, 2.607877, 2.280189, 2.204124, 2.260829,
};
const double ch21[] = {
6.053523, 3.777506, 3.562251, 3.601356, 3.569347, 3.410472,
};
const double ch22[] = {
11.932351, 4.803773, 3.782709, 3.390775, 3.226806,
3.218982, 3.382678, 3.773653, 4.788191, 11.905338,
};
const double ch23[] = {
7.145255, 5.637595, 4.049882, 3.025917, 2.326430, 1.773824, 1.235329,
};
const double ch24[] = {
2.268193, 2.372063, 2.323975, 2.034655, 1.736793,
};
const double ch25[] = {
2.231852, 2.383086, 2.341894, 2.031318, 1.729672, 1.486993,
};
const double ch26[] = {
3.993526, 2.308249, 2.075136, 2.038297, 2.036302, 2.078270, 2.295817, 4.017713,
};
const double ch27[] = {
2.205978,
2.175010,
2.215376,
2.473144,
};
const double ch28[] = {
2.321077, 2.371895, 2.338871, 2.057821, 1.755382,
};
const double ch29[] = {
2.222707, 2.372591, 2.301688, 1.991162, 1.695343,
};
const double ch30[] = {
2.599677,
2.026855,
2.138620,
2.229553,
};
const double ch31[] = {
5.791779, 2.636219, 2.103642, 1.967198, 1.962168, 2.096514, 2.641189, 5.780828,
};
const double ch32[] = {
2.006685, 2.545525, 2.477745, 2.335747, 2.194226, 2.076500,
};
applyCorrection(_hSL_topPt, ch15);
applyCorrection(_hSL_topPtTtbarSys, ch16);
applyCorrection(_hSL_topY, ch17);
applyCorrection(_hSL_ttbarDelPhi, ch18);
applyCorrection(_hSL_topPtLead, ch19);
applyCorrection(_hSL_topPtSubLead, ch20);
applyCorrection(_hSL_ttbarPt, ch21);
applyCorrection(_hSL_ttbarY, ch22);
applyCorrection(_hSL_ttbarMass, ch23);
applyCorrection(_hDL_topPt, ch24);
applyCorrection(_hDL_topPtTtbarSys, ch25);
applyCorrection(_hDL_topY, ch26);
applyCorrection(_hDL_ttbarDelPhi, ch27);
applyCorrection(_hDL_topPtLead, ch28);
applyCorrection(_hDL_topPtSubLead, ch29);
applyCorrection(_hDL_ttbarPt, ch30);
applyCorrection(_hDL_ttbarY, ch31);
applyCorrection(_hDL_ttbarMass, ch32);
}
if (_doShapeOnly) {
normalize(_hSL_topPt);
normalize(_hSL_topPtTtbarSys);
normalize(_hSL_topY);
normalize(_hSL_ttbarDelPhi);
normalize(_hSL_topPtLead);
normalize(_hSL_topPtSubLead);
normalize(_hSL_ttbarPt);
normalize(_hSL_ttbarY);
normalize(_hSL_ttbarMass);
normalize(_hDL_topPt);
normalize(_hDL_topPtTtbarSys);
normalize(_hDL_topY);
normalize(_hDL_ttbarDelPhi);
normalize(_hDL_topPtLead);
normalize(_hDL_topPtSubLead);
normalize(_hDL_ttbarPt);
normalize(_hDL_ttbarY);
normalize(_hDL_ttbarMass);
}
else {
const double s = 1. / sumOfWeights();
scale(_hSL_topPt, s);
scale(_hSL_topPtTtbarSys, s);
scale(_hSL_topY, s);
scale(_hSL_ttbarDelPhi, s);
scale(_hSL_topPtLead, s);
scale(_hSL_topPtSubLead, s);
scale(_hSL_ttbarPt, s);
scale(_hSL_ttbarY, s);
scale(_hSL_ttbarMass, s);
scale(_hDL_topPt, s);
scale(_hDL_topPtTtbarSys, s);
scale(_hDL_topY, s);
scale(_hDL_ttbarDelPhi, s);
scale(_hDL_topPtLead, s);
scale(_hDL_topPtSubLead, s);
scale(_hDL_ttbarPt, s);
scale(_hDL_ttbarY, s);
scale(_hDL_ttbarMass, s);
}
}
void applyCorrection(Histo1DPtr& h, const double* cf) {
for (auto& bin : h->bins()) {
bin.scaleW(cf[bin.index() - 1]);
}
}
private:
const bool _applyCorrection, _doShapeOnly;
Histo1DPtr _hSL_topPt, _hSL_topPtTtbarSys, _hSL_topY, _hSL_ttbarDelPhi, _hSL_topPtLead, _hSL_topPtSubLead,
_hSL_ttbarPt, _hSL_ttbarY, _hSL_ttbarMass;
Histo1DPtr _hDL_topPt, _hDL_topPtTtbarSys, _hDL_topY, _hDL_ttbarDelPhi, _hDL_topPtLead, _hDL_topPtSubLead,
_hDL_ttbarPt, _hDL_ttbarY, _hDL_ttbarMass;
};
RIVET_DECLARE_PLUGIN(CMS_2015_I1370682);
///////////////
// More PseudoTop implementation
namespace {
double PseudoTop::_tMass = 172.5 * GeV;
double PseudoTop::_wMass = 80.4 * GeV;
CmpState PseudoTop::compare(const Projection& p) const {
const PCmp fscmp = mkNamedPCmp(p, "FS");
if (fscmp != CmpState::EQ) return fscmp;
const PseudoTop& other = dynamic_cast<const PseudoTop&>(p);
CmpState cs_lepR = cmp(_lepR, other._lepR);
if (cs_lepR != CmpState::EQ) return cs_lepR;
CmpState cs_jetR = cmp(_jetR, other._jetR);
if (cs_jetR != CmpState::EQ) return cs_jetR;
CmpState cs_lepMinPt = cmp(_lepMinPt, other._lepMinPt);
if (cs_lepMinPt != CmpState::EQ) return cs_lepMinPt;
CmpState cs_jetMinPt = cmp(_jetMinPt, other._jetMinPt);
if (cs_jetMinPt != CmpState::EQ) return cs_jetMinPt;
CmpState cs_lepMaxEta = cmp(_lepMaxEta, other._lepMaxEta);
if (cs_lepMaxEta != CmpState::EQ) return cs_lepMaxEta;
CmpState cs_jetMaxEta = cmp(_jetMaxEta, other._jetMaxEta);
return cs_jetMaxEta;
}
void PseudoTop::cleanup(map<double, pair<size_t, size_t>>& v, const bool doCrossCleanup) const {
vector<map<double, pair<size_t, size_t>>::iterator> toErase;
set<size_t> usedLeg1, usedLeg2;
if (!doCrossCleanup) {
/// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) {
for (map<double, pair<size_t, size_t>>::iterator key = v.begin(); key != v.end(); ++key) {
const size_t leg1 = key->second.first;
const size_t leg2 = key->second.second;
if (usedLeg1.find(leg1) == usedLeg1.end() and usedLeg2.find(leg2) == usedLeg2.end()) {
usedLeg1.insert(leg1);
usedLeg2.insert(leg2);
}
else {
toErase.push_back(key);
}
}
}
else {
/// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) {
for (map<double, pair<size_t, size_t>>::iterator key = v.begin(); key != v.end(); ++key) {
const size_t leg1 = key->second.first;
const size_t leg2 = key->second.second;
if (usedLeg1.find(leg1) == usedLeg1.end() and usedLeg1.find(leg2) == usedLeg1.end()) {
usedLeg1.insert(leg1);
usedLeg1.insert(leg2);
}
else {
toErase.push_back(key);
}
}
}
/// @todo Reinstate when C++11 allowed: for (auto& key : toErase) v.erase(key);
for (size_t i = 0; i < toErase.size(); ++i) v.erase(toErase[i]);
}
void PseudoTop::project(const Event& e) {
// Leptons : do the lepton clustering anti-kt R=0.1 using stable photons and leptons not from hadron decay
// Neutrinos : neutrinos not from hadron decay
// MET : vector sum of all invisible particles in x-y plane
// Jets : anti-kt R=0.4 using all particles excluding neutrinos and particles used in lepton clustering
// add ghost B hadrons during the jet clustering to identify B jets.
// W->lv : dressed lepton and neutrino pairs
// W->jj : light flavored dijet
// W candidate : select lv or jj pairs which minimise |mW1-80.4|+|mW2-80.4|
// lepton-neutrino pair will be selected with higher priority
// t->Wb : W candidate + b jet
// t candidate : select Wb pairs which minimise |mtop1-172.5|+|mtop2-172.5|
_isValid = false;
_wDecays1.clear();
_wDecays2.clear();
_jets.clear();
_bjets.clear();
_ljets.clear();
_mode1 = _mode2 = CH_HADRON;
// Collect final state particles
Particles pForLep, pForJet;
Particles neutrinos; // Prompt neutrinos
/// @todo Avoid this unsafe jump into HepMC -- all this can be done properly via VisibleFS and HeavyHadrons projections
for (ConstGenParticlePtr p : HepMCUtils::particles(e.genEvent())) { //
const int status = p->status();
const int pid = p->pdg_id();
if (status == 1) {
Particle rp = *p;
if (!PID::isHadron(pid) && !rp.fromHadron()) {
// Collect particles not from hadron decay
if (rp.isNeutrino()) {
// Prompt neutrinos are kept in separate collection
neutrinos.push_back(rp);
}
else if (pid == PID::PHOTON || rp.isLepton()) {
// Leptons and photons for the dressing
pForLep.push_back(rp);
}
}
else if (!rp.isNeutrino()) {
// Use all particles from hadron decay
pForJet.push_back(rp);
}
}
else if (PID::isHadron(pid) && PID::hasBottom(pid)) {
// NOTE: Consider B hadrons with pT > 5GeV - not in CMS proposal
//if ( p->momentum().perp() < 5 ) continue;
// Do unstable particles, to be used in the ghost B clustering
// Use last B hadrons only
bool isLast = true;
for (ConstGenParticlePtr pp : HepMCUtils::particles(p->end_vertex(), Relatives::CHILDREN)) {
if (PID::hasBottom(pp->pdg_id())) {
isLast = false;
break;
}
}
if (!isLast) continue;
// Rescale momentum by 10^-20
/// @todo Why the factor of 1/rho() as well?
//Particle ghost(pdgId, FourMomentum(p->momentum())*1e-20/p->momentum().rho());
Particle ghost(pid, FourMomentum(p->momentum()));
ghost.setMomentum(ghost.momentum() * 1.e-20 / ghost.momentum().rho());
pForJet.push_back(ghost);
}
}
// Start object building from trivial thing - prompt neutrinos
isortByPt(neutrinos);
// Proceed to lepton dressing
FastJets fjLep(FinalState(), JetAlg::ANTIKT, _lepR);
fjLep.calc(pForLep);
Jets leptons;
vector<int> leptonsId;
set<int> dressedIdxs;
for (const Jet& lep : fjLep.jetsByPt(Cuts::pT > _lepMinPt && Cuts::abseta < _lepMaxEta)) {
double leadingPt = -1;
int leptonId = 0;
for (const Particle& p : lep.particles()) {
/// @warning Barcodes aren't future-proof in HepMC
dressedIdxs.insert(HepMCUtils::uniqueId(p.genParticle()));
if (p.isLepton() && p.pT() > leadingPt) {
leadingPt = p.pT();
leptonId = p.pid();
}
}
if (leptonId == 0) continue;
leptons.push_back(lep);
leptonsId.push_back(leptonId);
}
// Re-use particles not used in lepton dressing
for (const Particle& rp : pForLep) {
const int barcode = HepMCUtils::uniqueId(rp.genParticle());
// Skip if the particle is used in dressing
if (dressedIdxs.find(barcode) != dressedIdxs.end()) continue;
// Put back to be used in jet clustering
pForJet.push_back(rp);
}
// Then do the jet clustering
FastJets fjJet(FinalState(), JetAlg::ANTIKT, _jetR);
//fjJet.useInvisibles(); // NOTE: CMS proposal to remove neutrinos (AB: wouldn't work anyway, since they were excluded from clustering inputs)
fjJet.calc(pForJet);
for (const Jet& jet : fjJet.jetsByPt(Cuts::pT > _jetMinPt && Cuts::abseta < _jetMaxEta)) {
_jets.push_back(jet);
bool isBJet = false;
for (const Particle& rp : jet.particles()) {
if (PID::hasBottom(rp.pid())) {
isBJet = true;
break;
}
}
if (isBJet)
_bjets.push_back(jet);
else
_ljets.push_back(jet);
}
// Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
if (_bjets.size() < 2) return; // Ignore single top for now
map<double, pair<size_t, size_t>> wLepCandIdxs;
map<double, pair<size_t, size_t>> wHadCandIdxs;
// Collect leptonic-decaying W's
for (size_t iLep = 0, nLep = leptons.size(); iLep < nLep; ++iLep) {
const Jet& lep = leptons.at(iLep);
for (size_t iNu = 0, nNu = neutrinos.size(); iNu < nNu; ++iNu) {
const Particle& nu = neutrinos.at(iNu);
const double m = (lep.momentum() + nu.momentum()).mass();
const double dm = std::abs(m - _wMass);
wLepCandIdxs[dm] = make_pair(iLep, iNu);
}
}
// Continue to hadronic decaying W's
for (size_t i = 0, nLjet = _ljets.size(); i < nLjet; ++i) {
const Jet& ljet1 = _ljets[i];
for (size_t j = i + 1; j < nLjet; ++j) {
const Jet& ljet2 = _ljets[j];
const double m = (ljet1.momentum() + ljet2.momentum()).mass();
const double dm = std::abs(m - _wMass);
wHadCandIdxs[dm] = make_pair(i, j);
}
}
// Cleanup W candidate, choose pairs with minimum dm if they share decay products
cleanup(wLepCandIdxs);
cleanup(wHadCandIdxs, true);
const size_t nWLepCand = wLepCandIdxs.size();
const size_t nWHadCand = wHadCandIdxs.size();
if (nWLepCand + nWHadCand < 2) return; // We skip single top
int w1Q = 1, w2Q = -1;
int w1dau1Id = 1, w2dau1Id = -1;
FourMomentum w1dau1LVec, w1dau2LVec;
FourMomentum w2dau1LVec, w2dau2LVec;
if (nWLepCand == 0) { // Full hadronic case
const pair<size_t, size_t>& idPair1 = wHadCandIdxs.begin()->second;
const pair<size_t, size_t>& idPair2 = (++wHadCandIdxs.begin())->second; ///< @todo Reinstate std::next
const Jet& w1dau1 = _ljets[idPair1.first];
const Jet& w1dau2 = _ljets[idPair1.second];
const Jet& w2dau1 = _ljets[idPair2.first];
const Jet& w2dau2 = _ljets[idPair2.second];
w1dau1LVec = w1dau1.momentum();
w1dau2LVec = w1dau2.momentum();
w2dau1LVec = w2dau1.momentum();
w2dau2LVec = w2dau2.momentum();
}
else if (nWLepCand == 1) { // Semi-leptonic case
const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
const pair<size_t, size_t>& idPair2 = wHadCandIdxs.begin()->second;
const Jet& w1dau1 = leptons[idPair1.first];
const Particle& w1dau2 = neutrinos[idPair1.second];
const Jet& w2dau1 = _ljets[idPair2.first];
const Jet& w2dau2 = _ljets[idPair2.second];
w1dau1LVec = w1dau1.momentum();
w1dau2LVec = w1dau2.momentum();
w2dau1LVec = w2dau1.momentum();
w2dau2LVec = w2dau2.momentum();
w1dau1Id = leptonsId[idPair1.first];
w1Q = w1dau1Id > 0 ? -1 : 1;
w2Q = -w1Q;
switch (w1dau1Id) {
case 13:
case -13: _mode1 = CH_MUON; break;
case 11:
case -11: _mode1 = CH_ELECTRON; break;
}
}
else { // Full leptonic case
const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
const pair<size_t, size_t>& idPair2 = (++wLepCandIdxs.begin())->second; ///< @todo Reinstate std::next
const Jet& w1dau1 = leptons[idPair1.first];
const Particle& w1dau2 = neutrinos[idPair1.second];
const Jet& w2dau1 = leptons[idPair2.first];
const Particle& w2dau2 = neutrinos[idPair2.second];
w1dau1LVec = w1dau1.momentum();
w1dau2LVec = w1dau2.momentum();
w2dau1LVec = w2dau1.momentum();
w2dau2LVec = w2dau2.momentum();
w1dau1Id = leptonsId[idPair1.first];
w2dau1Id = leptonsId[idPair2.first];
w1Q = w1dau1Id > 0 ? -1 : 1;
w2Q = w2dau1Id > 0 ? -1 : 1;
switch (w1dau1Id) {
case 13:
case -13: _mode1 = CH_MUON; break;
case 11:
case -11: _mode1 = CH_ELECTRON; break;
}
switch (w2dau1Id) {
case 13:
case -13: _mode2 = CH_MUON; break;
case 11:
case -11: _mode2 = CH_ELECTRON; break;
}
}
const FourMomentum w1LVec = w1dau1LVec + w1dau2LVec;
const FourMomentum w2LVec = w2dau1LVec + w2dau2LVec;
// Combine b jets
double sumDm = 1e9;
FourMomentum b1LVec, b2LVec;
for (size_t i = 0, n = _bjets.size(); i < n; ++i) {
const Jet& bjet1 = _bjets[i];
const double mtop1 = (w1LVec + bjet1.momentum()).mass();
const double dmtop1 = std::abs(mtop1 - _tMass);
for (size_t j = 0; j < n; ++j) {
if (i == j) continue;
const Jet& bjet2 = _bjets[j];
const double mtop2 = (w2LVec + bjet2.momentum()).mass();
const double dmtop2 = std::abs(mtop2 - _tMass);
if (sumDm <= dmtop1 + dmtop2) continue;
sumDm = dmtop1 + dmtop2;
b1LVec = bjet1.momentum();
b2LVec = bjet2.momentum();
}
}
if (sumDm >= 1e9) return; // Failed to make top, but this should not happen.
const FourMomentum t1LVec = w1LVec + b1LVec;
const FourMomentum t2LVec = w2LVec + b2LVec;
// Put all of them into candidate collection
_t1 = Particle(w1Q * 6, t1LVec);
_b1 = Particle(w1Q * 5, b1LVec);
_wDecays1.push_back(Particle(w1dau1Id, w1dau1LVec));
_wDecays1.push_back(Particle(-w1dau1Id + w1Q, w1dau2LVec));
_t2 = Particle(w2Q * 6, t2LVec);
_b2 = Particle(w2Q * 5, b2LVec);
_wDecays2.push_back(Particle(w2dau1Id, w2dau1LVec));
_wDecays2.push_back(Particle(-w2dau1Id + w2Q, w2dau2LVec));
_isValid = true;
}
}
}