Rivet analyses referenceATLAS_2017_I1614149Resolved and boosted ttbar l+jets cross sections at 13 TeVExperiment: ATLAS (LHC) Inspire ID: 1614149 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Measurements of differential cross-sections of top-quark pair production in fiducial phase-spaces are presented as a function of top-quark and $t\bar{t}$ system kinematic observables in proton-proton collisions at a centre-of-mass energy of $\sqrt{s} = 13$ TeV. The data set corresponds to an integrated luminosity of 3.2 fb${}^{-1}$, recorded in 2015 with the ATLAS detector at the CERN Large Hadron Collider. Events with exactly one electron or muon and at least two jets in the final state are used for the measurement. Two separate selections are applied that each focus on different top-quark momentum regions, referred to as resolved and boosted topologies of the $t\bar{t}$ final state. The measured spectra are corrected for detector effects and are compared to several Monte Carlo simulations by means of calculated $\chi^2$ and $p$-values. Source code: ATLAS_2017_I1614149.cc 1// -*- C++ -*
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/VetoedFinalState.hh"
5#include "Rivet/Projections/IdentifiedFinalState.hh"
6#include "Rivet/Projections/PromptFinalState.hh"
7#include "Rivet/Projections/LeptonFinder.hh"
8#include "Rivet/Projections/FastJets.hh"
9#include "Rivet/Projections/MissingMomentum.hh"
10
11#include "fastjet/contrib/Njettiness.hh"
12#include "fastjet/contrib/Nsubjettiness.hh"
13#include "fastjet/contrib/NjettinessPlugin.hh"
14
15namespace Rivet {
16
17
18 class ATLAS_2017_I1614149 : public Analysis {
19 public:
20
21 /// Constructor
22 ///@brief: Resolved and boosted ttbar l+jets cross sections at 13 TeV
23 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1614149);
24
25 void init() {
26 // Eta ranges
27 Cut eta_full = (Cuts::abseta < 5.0);
28 Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV);
29
30 // All final state particles
31 FinalState fs(eta_full);
32
33 IdentifiedFinalState all_photons(fs);
34 all_photons.acceptIdPair(PID::PHOTON);
35
36 // Get photons to dress leptons
37 IdentifiedFinalState ph_id(fs);
38 ph_id.acceptIdPair(PID::PHOTON);
39
40 // Projection to find the electrons
41 IdentifiedFinalState el_id(fs);
42 el_id.acceptIdPair(PID::ELECTRON);
43
44 PromptFinalState photons(ph_id);
45 photons.acceptTauDecays(true);
46 declare(photons, "photons");
47
48 PromptFinalState electrons(el_id);
49 electrons.acceptTauDecays(true);
50 LeptonFinder dressedelectrons(electrons, photons, 0.1, lep_cuts);
51 declare(dressedelectrons, "elecs");
52 LeptonFinder ewdressedelectrons(electrons, all_photons, 0.1, eta_full);
53
54 // Projection to find the muons
55 IdentifiedFinalState mu_id(fs);
56 mu_id.acceptIdPair(PID::MUON);
57
58 PromptFinalState muons(mu_id);
59 muons.acceptTauDecays(true);
60 LeptonFinder dressedmuons(muons, photons, 0.1, lep_cuts);
61 declare(dressedmuons, "muons");
62 LeptonFinder ewdressedmuons(muons, all_photons, 0.1, eta_full);
63
64 // Projection to find MET
65 declare(MissingMomentum(fs), "MET");
66
67 // remove prompt neutrinos from jet clustering
68 IdentifiedFinalState nu_id(fs);
69 nu_id.acceptNeutrinos();
70 PromptFinalState neutrinos(nu_id);
71 neutrinos.acceptTauDecays(true);
72
73 // Jet clustering.
74 VetoedFinalState vfs(fs);
75 vfs.addVetoOnThisFinalState(ewdressedelectrons);
76 vfs.addVetoOnThisFinalState(ewdressedmuons);
77 vfs.addVetoOnThisFinalState(neutrinos);
78 FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL);
79 declare(jets, "jets");
80
81 // Addition of the large-R jets
82 VetoedFinalState vfs1(fs);
83 vfs1.addVetoOnThisFinalState(neutrinos);
84 FastJets fjets(vfs1, JetAlg::ANTIKT, 1.);
85 fjets.useInvisibles(JetInvisibles::NONE);
86 fjets.useMuons(JetMuons::NONE);
87 declare(fjets, "fjets");
88
89 bookHists("top_pt_res", 15);
90 bookHists("top_absrap_res", 17);
91 bookHists("ttbar_pt_res", 19);
92 bookHists("ttbar_absrap_res", 21);
93 bookHists("ttbar_m_res", 23);
94 bookHists("top_pt_boost", 25);
95 bookHists("top_absrap_boost", 27);
96
97 }
98
99
100 void analyze(const Event& event) {
101
102 // Get the selected objects, using the projections.
103 DressedLeptons electrons = apply<LeptonFinder>(event, "elecs").dressedLeptons();
104 DressedLeptons muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
105 const Jets& jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
106 const PseudoJets& all_fjets = apply<FastJets>(event, "fjets").pseudojetsByPt();
107
108 // get MET
109 const Vector3 met = apply<MissingMomentum>(event, "MET").vectorMPT();
110
111 Jets bjets, lightjets;
112 for (const Jet& jet : jets) {
113 bool b_tagged = jet.bTags(Cuts::pT > 5*GeV).size();
114 if ( b_tagged && bjets.size() < 2) bjets +=jet;
115 else lightjets += jet;
116 }
117
118 // Implementing large-R jets definition
119 // trim the jets
120 PseudoJets trimmed_fatJets;
121 float Rfilt = 0.2;
122 float pt_fraction_min = 0.05;
123 fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, Rfilt), fastjet::SelectorPtFractionMin(pt_fraction_min));
124 for (PseudoJet pjet : all_fjets) trimmed_fatJets += trimmer(pjet);
125 trimmed_fatJets = fastjet::sorted_by_pt(trimmed_fatJets);
126 PseudoJets trimmed_jets;
127 for (unsigned int i = 0; i < trimmed_fatJets.size(); ++i) {
128 FourMomentum tj_mom = momentum(trimmed_fatJets[i]);
129 if (tj_mom.pt() <= 300*GeV) continue;
130 if (tj_mom.abseta() >= 2.0) continue;
131 trimmed_jets.push_back(trimmed_fatJets[i]);
132 }
133
134 bool single_electron = (electrons.size() == 1) && (muons.empty());
135 bool single_muon = (muons.size() == 1) && (electrons.empty());
136
137 DressedLepton *lepton = NULL;
138 if (single_electron) lepton = &electrons[0];
139 else if (single_muon) lepton = &muons[0];
140
141 if (!single_electron && !single_muon) vetoEvent;
142
143 bool pass_resolved = true;
144 bool num_b_tagged_jets = (bjets.size() == 2);
145 if (!num_b_tagged_jets) pass_resolved = false;
146
147 if (jets.size() < 4) pass_resolved = false;
148
149 bool pass_boosted = true;
150 int fatJetIndex = -1;
151 bool passTopTag = false;
152 bool passDphi = false;
153 bool passAddJet = false;
154 bool goodLepJet = false;
155 bool lepbtag = false;
156 bool hadbtag=false;
157 vector<int> lepJetIndex;
158 vector<int> jet_farFromHadTopJetCandidate;
159 if (met.mod() < 20*GeV) pass_boosted = false;
160 if (pass_boosted) {
161 double transmass = _mT(lepton->momentum(), met);
162 if (transmass + met.mod() < 60*GeV) pass_boosted = false;
163 }
164 if (pass_boosted) {
165 if (trimmed_jets.size() >= 1) {
166 for (unsigned int j = 0; j<trimmed_jets.size(); ++j) {
167 if (tau32( trimmed_jets.at(j), 1. ) < 0.75 &&
168 momentum(trimmed_jets.at(j)).mass() > 100*GeV &&
169 momentum(trimmed_jets.at(j)).pt() > 300*GeV &&
170 momentum(trimmed_jets.at(j)).pt() < 1500*GeV &&
171 fabs(momentum(trimmed_jets.at(j)).eta()) < 2.) {
172 passTopTag = true;
173 fatJetIndex = j;
174 break;
175 }
176 }
177 }
178 }
179 if(!passTopTag && fatJetIndex == -1) pass_boosted = false;
180 if (pass_boosted) {
181 double dPhi_fatjet = deltaPhi(lepton->phi(), momentum(trimmed_jets.at(fatJetIndex)).phi());
182 double dPhi_fatjet_lep_cut = 1.0; //2.3
183 if (dPhi_fatjet > dPhi_fatjet_lep_cut ) {
184 passDphi = true;
185 }
186 }
187 if (!passDphi) pass_boosted = false;
188 if (bjets.empty()) pass_boosted = false;
189 if (pass_boosted) {
190 for (unsigned int sj = 0; sj < jets.size(); ++sj) {
191 double dR = deltaR(jets.at(sj).momentum(), momentum(trimmed_jets.at(fatJetIndex)));
192 if(dR > 1.5) {
193 passAddJet = true;
194 jet_farFromHadTopJetCandidate.push_back(sj);
195 }
196 }
197 }
198 if (!passAddJet) pass_boosted = false;
199 if (pass_boosted) {
200 for (int ltj : jet_farFromHadTopJetCandidate) {
201 double dR_jet_lep = deltaR(jets.at(ltj).momentum(), lepton->momentum());
202 double dR_jet_lep_cut = 2.0;//1.5
203 if (dR_jet_lep < dR_jet_lep_cut) {
204 lepJetIndex.push_back(ltj);
205 goodLepJet = true;
206 }
207 }
208 }
209 if(!goodLepJet) pass_boosted = false;
210 if (pass_boosted) {
211 for (int lepj : lepJetIndex) {
212 lepbtag = jets.at(lepj).bTags(Cuts::pT > 5*GeV).size();
213 if (lepbtag) break;
214 }
215 }
216 double dR_fatBjet_cut = 1.0;
217 if (pass_boosted) {
218 for (const Jet& bjet : bjets) {
219 hadbtag |= deltaR(momentum(trimmed_jets.at(fatJetIndex)), bjet) < dR_fatBjet_cut;
220 }
221 }
222
223 if (!(lepbtag || hadbtag)) pass_boosted = false;
224
225 FourMomentum pbjet1; //Momentum of bjet1
226 FourMomentum pbjet2; //Momentum of bjet
227 int Wj1index = -1, Wj2index = -1;
228
229 if (pass_resolved) {
230
231 if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) {
232 pbjet1 = bjets[0].momentum();
233 pbjet2 = bjets[1].momentum();
234 } else {
235 pbjet1 = bjets[1].momentum();
236 pbjet2 = bjets[0].momentum();
237 }
238
239 double bestWmass = 1000.0*TeV;
240 double mWPDG = 80.399*GeV;
241 for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) {
242 for (unsigned int j = i + 1; j < lightjets.size(); ++j) {
243 double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass();
244 if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) {
245 bestWmass = wmass;
246 Wj1index = i;
247 Wj2index = j;
248 }
249 }
250 }
251
252 FourMomentum pjet1 = lightjets[Wj1index].momentum();
253 FourMomentum pjet2 = lightjets[Wj2index].momentum();
254
255 // compute hadronic W boson
256 FourMomentum pWhadron = pjet1 + pjet2;
257 double pz = computeneutrinoz(lepton->momentum(), met);
258 FourMomentum ppseudoneutrino( sqrt(sqr(met.x()) + sqr(met.y()) + sqr(pz)), met.x(), met.y(), pz);
259
260 //compute leptonic, hadronic, combined pseudo-top
261 FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1;
262 FourMomentum ppseudotophadron = pbjet2 + pWhadron;
263 FourMomentum pttbar = ppseudotoplepton + ppseudotophadron;
264
265 fillHists("top_pt_res", ppseudotophadron.pt()/GeV);
266 fillHists("top_absrap_res", ppseudotophadron.absrap());
267 fillHists("ttbar_pt_res", pttbar.pt()/GeV);
268 fillHists("ttbar_absrap_res", pttbar.absrap());
269 fillHists("ttbar_m_res", pttbar.mass()/GeV);
270 }
271
272 if (pass_boosted) {// Boosted selection
273 double hadtop_pt= momentum(trimmed_jets.at(fatJetIndex)).pt() / GeV;
274 double hadtop_absrap= momentum(trimmed_jets.at(fatJetIndex)).absrap();
275 fillHists("top_pt_boost", hadtop_pt);
276 fillHists("top_absrap_boost", hadtop_absrap);
277 }
278 }
279
280
281 void finalize() {
282 // Normalize to cross-section
283 const double sf = (crossSection()/picobarn / sumOfWeights());
284 for (HistoMap::value_type& hist : _h) {
285 scale(hist.second, sf);
286 if (hist.first.find("_norm") != string::npos) normalize(hist.second);
287 }
288 }
289
290
291 void bookHists(std::string name, unsigned int index) {
292 book(_h[name], index, 1 ,1);
293 book(_h[name + "_norm"], index + 1, 1, 1);
294 }
295
296
297 void fillHists(std::string name, double value) {
298 _h[name]->fill(value);
299 _h[name + "_norm"]->fill(value);
300 }
301
302
303 double _mT(const FourMomentum &l, const Vector3 &met) const {
304 return sqrt(2.0 * l.pT() * met.mod() * (1 - cos(deltaPhi(l, met))) );
305 }
306
307
308 double tau32(const fastjet::PseudoJet &jet, double jet_rad) const {
309 double alpha = 1.0;
310 fjcontrib::NormalizedCutoffMeasure normalized_measure(alpha, jet_rad, 1000000);
311 // WTA definition
312 // Nsubjettiness::OnePass_WTA_KT_Axes wta_kt_axes;
313 // as in JetSubStructure recommendations
314 fjcontrib::KT_Axes kt_axes;
315
316 /// NsubjettinessRatio uses the results from Nsubjettiness to calculate the ratio
317 /// tau_N/tau_M, where N and M are specified by the user. The ratio of different tau values
318 /// is often used in analyses, so this class is helpful to streamline code.
319 fjcontrib::NsubjettinessRatio tau32_kt(3, 2, kt_axes, normalized_measure);
320
321 double tau32 = tau32_kt.result(jet);
322 return tau32;
323 }
324
325
326 double computeneutrinoz(const FourMomentum& lepton, const Vector3 &met) const {
327 //computing z component of neutrino momentum given lepton and met
328 double pzneutrino;
329 double m_W = 80.399; // in GeV, given in the paper
330 double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.x() + lepton.py() * met.y());
331 double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
332 double b = -2*k*lepton.pz();
333 double c = sqr( lepton.E() ) * sqr( met.mod() ) - sqr( k );
334 double discriminant = sqr(b) - 4 * a * c;
335 double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
336 if (discriminant < 0) pzneutrino = - b / (2 * a); //if the discriminant is negative
337 else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
338 double absquad[2];
339 for (int n=0; n<2; ++n) absquad[n] = fabs(quad[n]);
340 if (absquad[0] < absquad[1]) pzneutrino = quad[0];
341 else pzneutrino = quad[1];
342 }
343
344 return pzneutrino;
345 }
346
347
348 private:
349
350 /// @name Objects that are used by the event selection decisions
351 typedef map<string, Histo1DPtr> HistoMap;
352 HistoMap _h;
353
354 };
355
356
357 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1614149);
358
359}
|