Rivet analyses referenceATLAS_2021_I1941095energy asymmetry in ttj @ 13 TeVExperiment: ATLAS (LHC) Inspire ID: 1941095 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
A measurement of the energy asymmetry in jet-associated top-quark pair production is presented using 139 fb$^{-1}$ of data collected by the ATLAS detector at the Large Hadron Collider during $pp$ collisions at $\sqrt{s}=13$ TeV. The observable measures the different probability of top and antitop quarks to have the higher energy as a function of the jet scattering angle with respect to the beam axis. The energy asymmetry is measured in the semileptonic $t\bar{t}$ decay channel, and the hadronically decaying top quark must have transverse momentum above 350 GeV. The results are corrected for detector effects to particle level in three bins of the scattering angle of the associated jet. The measurement agrees with the SM prediction at next-to-leading-order accuracy in quantum chromodynamics in all three bins. In the bin with the largest expected asymmetry, where the jet is emitted perpendicular to the beam, the energy asymmetry is measured to be $-0.043\pm0.020$, in agreement with the SM prediction of $-0.037\pm0.003$. Interpreting this result in the framework of the Standard Model effective field theory (SMEFT), it is shown that the energy asymmetry is sensitive to the top-quark chirality in four-quark operators and is therefore a valuable new observable in global SMEFT fits. Source code: ATLAS_2021_I1941095.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/InvisibleFinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/LeptonFinder.hh"
7#include "Rivet/Projections/MissingMomentum.hh"
8#include "Rivet/Projections/PromptFinalState.hh"
9#include "Rivet/Projections/PartonicTops.hh"
10#include "Rivet/Projections/VetoedFinalState.hh"
11#include "Rivet/Tools/MendelMin.hh"
12#include "fastjet/tools/Filter.hh"
13
14namespace Rivet {
15
16
17 /// @brief Energy asymmetry in ttj at 13 TeV
18 class ATLAS_2021_I1941095 : public Analysis {
19 public:
20
21 /// Constructor
22 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2021_I1941095);
23
24
25 /// @name Analysis methods
26 /// @{
27
28 /// Book histograms and initialise projections before the run
29 void init() {
30
31 // Declare projections
32
33 // Photons
34 PromptFinalState promptphotons(Cuts::abspid == PID::PHOTON, TauDecaysAs::NONPROMPT);
35
36 // Electrons
37 PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
38 LeptonFinder all_dressed_el(bare_el, promptphotons, 0.1, Cuts::abseta < 2.5);
39 LeptonFinder electrons(bare_el, promptphotons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
40 declare(electrons,"electrons");
41
42 // Muons
43 PromptFinalState bare_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
44 LeptonFinder all_dressed_mu(bare_mu, promptphotons, 0.1, Cuts::abseta < 2.5);
45 LeptonFinder muons(bare_mu, promptphotons, 0.1, Cuts::abseta <2.5 && Cuts::pT > 25*GeV);
46 declare(muons,"muons");
47
48 // AntiKt4TruthWZJets as AntiKt4TruthWZJets, but w/o photons from hadrons in dressing
49 const InvisibleFinalState invisibles(OnlyPrompt::YES, TauDecaysAs::PROMPT);
50 VetoedFinalState vfs(FinalState(Cuts::abseta < 5.0)); // changed from 4.5 to 5.0
51 vfs.addVetoOnThisFinalState(all_dressed_el);
52 vfs.addVetoOnThisFinalState(all_dressed_mu);
53 vfs.addVetoOnThisFinalState(invisibles); // new
54 FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL); // changed invisible from DECAY to ALL
55 declare(jets,"jets");
56
57 // AntiKt10TruthTrimmedPtFrac5SmallR20Jets
58 FinalState fs(Cuts::abseta < 5.0);
59 FastJets fjets(fs, JetAlg::ANTIKT, 1.0, JetMuons::NONE, JetInvisibles::NONE);
60 _trimmer = fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2),
61 fastjet::SelectorPtFractionMin(0.05));
62 declare(fjets,"fjets");
63
64 // Missing momentum
65 declare(MissingMomentum(), "MissingMomentum");
66
67 // Parton level top quarks after FSR
68 // options are: decaymode, emu_from_prompt_tau, include_hadronic_taus
69 declare(PartonicTops(TopDecay::E_MU, PromptEMuFromTau::YES, InclHadronicTau::NO), "PartonicTops_EMU");
70 declare(PartonicTops(TopDecay::E_MU, PromptEMuFromTau::NO, InclHadronicTau::NO), "PartonicTops_EMU_notau");
71 declare(PartonicTops(TopDecay::HADRONIC, PromptEMuFromTau::NO, InclHadronicTau::YES), "PartonicTops_HADRONIC");
72 declare(PartonicTops(TopDecay::HADRONIC, PromptEMuFromTau::NO, InclHadronicTau::NO), "PartonicTops_HADRONIC_notau");
73
74 // Book histograms
75 const Estimate1D& ref_asymm = refData(1, 1, 1);
76 book(_h["pos"], "_thetaj_opt_depos", ref_asymm.xEdges());
77 book(_h["neg"], "_thetaj_opt_deneg", ref_asymm.xEdges());
78 book(_asymm, 1, 1, 1);
79
80 }
81
82
83 /// Perform the per-event analysis
84 void analyze(const Event& event) {
85
86 // Parton-level top quarks // after FSR
87 const Particles partonicTops_EMU = apply<ParticleFinder>(event, "PartonicTops_EMU").particlesByPt();
88 const Particles partonicTops_EMU_notau = apply<ParticleFinder>(event, "PartonicTops_EMU_notau").particlesByPt();
89 const Particles partonicTops_HADRONIC = apply<ParticleFinder>(event, "PartonicTops_HADRONIC").particlesByPt();
90 const Particles partonicTops_HADRONIC_notau = apply<ParticleFinder>(event, "PartonicTops_HADRONIC_notau").particlesByPt();
91
92 // Filter semi-leptonic (e,mu,tau) events: Veto dileptonic/nonleptonic events and events with 2 taus
93 int nLeptons = partonicTops_EMU.size() + partonicTops_HADRONIC.size() - partonicTops_HADRONIC_notau.size();
94 if (nLeptons != 1) vetoEvent;
95
96 // Get the selected objects, using the projections.
97 DressedLeptons electrons = apply<LeptonFinder>(event, "electrons").dressedLeptons();
98 DressedLeptons muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
99 const Jets& jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
100 const Jets& fjets = apply<FastJets>(event, "fjets").jetsByPt(Cuts::pT > 200*GeV && Cuts::abseta < 2.0);
101 PseudoJets ljets;
102 for (const Jet& fjet : fjets) { ljets += _trimmer(fjet); }
103 sort(ljets.begin(), ljets.end(), [](PseudoJet const &l, PseudoJet const &r) { return l.pt() > r.pt(); });
104 const FourMomentum& met = apply<MissingMomentum>(event, "MissingMomentum").missingMomentum();
105
106 // Overlap removal
107 for (const Jet& jet : jets) {
108 idiscard(electrons, deltaRLess(jet, 0.4, RAPIDITY));
109 idiscard(muons, deltaRLess(jet, 0.4, RAPIDITY));
110 }
111
112 // Reconstruct event
113
114 // Lepton l
115 size_t n_el_25 = 0, n_el_27 = 0;
116 for (const DressedLepton& electron : electrons ) {
117 if (electron.pT() >= 25*GeV) ++n_el_25;
118 if (electron.pT() >= 27*GeV) ++n_el_27;
119 }
120 size_t n_mu_25 = 0, n_mu_27 = 0;
121 for (const DressedLepton& muon : muons ) {
122 if (muon.pT() >= 25*GeV) ++n_mu_25;
123 if (muon.pT() >= 27*GeV) ++n_mu_27;
124 }
125 if ((n_el_25 + n_mu_25 != 1) || (n_el_27 + n_mu_27 != 1)) vetoEvent;
126 DressedLepton lepton = (n_el_27 == 1)? electrons[0] : muons[0];
127 FourMomentum l = lepton.mom();
128 int lep_charge = lepton.charge();
129
130 // Neutrino nu
131 FourMomentum nu = getNeutrino(l, met);
132
133 // Hadronic top candidate jh
134 int jh_idx = -1;
135 for (size_t ijet = 0; ijet < ljets.size(); ++ijet) {
136 Jet ljet = Jet(ljets[ijet]);
137 if (ljet.pT() < 350*GeV) continue;
138 if (ljet.abseta() > 2.0) continue;
139 if (ljet.mass() < 140*GeV) continue;
140 if (deltaPhi(ljet,l) < 1.0) continue;
141 bool btagged = false;
142 for (const Jet& jet : jets) {
143 if ( jet.bTagged(Cuts::pT > 5*GeV) ) {
144 if ( deltaR(ljet, jet) < 1.0 ) {
145 btagged = true;
146 break;
147 }
148 }
149 }
150 if (btagged) {
151 jh_idx = ijet;
152 break;
153 }
154 }
155 if ( jh_idx == -1 ) vetoEvent;
156 FourMomentum jh = Jet(ljets[jh_idx]).mom();
157
158 // Leptonic top b-jet candidate jl
159 int jl_idx = -1;
160 for (size_t ijet = 0; ijet < jets.size(); ++ijet) {
161 if ( !jets[ijet].bTagged(Cuts::pT > 5*GeV) ) continue;
162 if ( deltaR(jets[ijet], l) > 2.0 ) continue;
163 if ( deltaR(jets[ijet], jh) < 1.5 ) continue;
164 jl_idx = ijet;
165 break;
166 }
167 if ( jl_idx == -1) {
168 for (size_t ijet = 0; ijet < jets.size(); ++ijet) {
169 if ( deltaR(jets[ijet], l) > 2.0 ) continue;
170 if ( deltaR(jets[ijet], jh) < 1.5 ) continue;
171 jl_idx = ijet;
172 break;
173 }
174 }
175 if ( jl_idx == -1 ) vetoEvent;
176 FourMomentum jl = Jet(jets[jl_idx]).mom();
177
178 // b-tagging
179 size_t n_btagged = 0;
180 size_t n_btagged_matched = 1; // Large-jet jh is b-tagged
181 for (const Jet& jet : jets ) {
182 if (jet.bTagged(Cuts::pT > 5*GeV)) ++n_btagged;
183 }
184 if ( jets[jl_idx].bTagged(Cuts::pT > 5*GeV) ) ++n_btagged_matched;
185 if ( n_btagged >= 2 && n_btagged_matched < 2 ) vetoEvent;
186
187 // Associated jet candidate ja
188 int ja_idx = -1;
189 for (int ijet = 0; ijet < int(jets.size()); ++ijet) {
190 if ( ijet == jl_idx ) continue;
191 if ( jets[ijet].pT() < 100*GeV ) continue;
192 if ( deltaR(jets[ijet], jh) < 1.5 ) continue;
193 if ( deltaR(jets[ijet], l) < 0.4 ) continue;
194 ja_idx = ijet;
195 break;
196 }
197 if ( ja_idx == -1 ) vetoEvent;
198 FourMomentum ja = jets[ja_idx].mom();
199
200 FourMomentum thad = jh;
201 FourMomentum tlep = l+nu+jl;
202 FourMomentum top = lep_charge > 0 ? tlep : thad;
203 FourMomentum tbar = lep_charge > 0 ? thad : tlep;
204 FourMomentum ttbar = top + tbar;
205 FourMomentum ttj = top + tbar + ja;
206
207 // Boost into ttj reference frame
208 FourMomentum ttj_inv( ttj.E(), -ttj.px(), -ttj.py(), -ttj.pz() );
209 Vector3 boostVector = ttj_inv.betaVec();
210 LorentzTransform lt_boost;
211 lt_boost.setBetaVec( boostVector );
212 FourMomentum top_boosted = lt_boost.transform( top );
213 FourMomentum tbar_boosted = lt_boost.transform( tbar );
214 FourMomentum ja_boosted = lt_boost.transform( ja );
215
216 // Get observables
217 const double deltaE = top_boosted.E() - tbar_boosted.E();
218 const double thetaj_opt = ttj.rapidity() > 0 ? ja_boosted.theta() : pi - ja_boosted.theta();
219
220 // Fill auxiliary histograms
221 _h[deltaE > 0 ? "pos" : "neg"]->fill(thetaj_opt);
222
223 }
224
225
226 /// Normalise histograms etc., after the run
227 void finalize() {
228
229 scale(_h, crossSection()/picobarn/sumW());
230
231 // Calculate differential energy asymmetry
232 asymm(_h["pos"], _h["neg"], _asymm);
233
234 }
235
236 /// @}
237
238
239 private:
240
241 fastjet::Filter _trimmer;
242
243 // Histograms
244 map<string, Histo1DPtr> _h;
245 Estimate1DPtr _asymm;
246
247
248 static double delta2_fcn(const MendelMin::Params& p, const MendelMin::Params& pfix) {
249 double delta2 = 0;
250 double alpha = p[0]*6.30-3.15; // Map p[0] in [0,1] to alpha in [-3.15,3.15]
251 double r = pfix[0];
252 double dphi = pfix[1];
253 double l_pt = pfix[2];
254 double l_m = pfix[3];
255 double n_px = pfix[4];
256 double n_py = pfix[5];
257 r /= sqrt(l_pt * l_pt + l_m * l_m) - l_pt * cos(dphi + alpha);
258 FourMomentum neut(0.0, n_px, n_py, 0.0); // E, px, py, pz
259 neut.setE(neut.p());
260 FourMomentum neut_new( 0.0, r * neut.p() * cos(neut.phi() + alpha), r * neut.p() * sin(neut.phi() + alpha), 0.0 );
261 neut_new.setE(neut_new.p());
262 delta2 = pow((neut_new.px() - neut.px()), 2) + pow((neut_new.py() - neut.py()), 2);
263 return delta2;
264 }
265
266
267 FourMomentum getNeutrino(const FourMomentum& lepton, const FourMomentum& met) {
268 const double m_mWpdg = 80.4*GeV;
269 double pxNu = met.px();
270 double pyNu = met.py();
271 double ptNu = met.pt();
272 double pzNu;
273
274 double c1 = pow(m_mWpdg,2) - pow(lepton.mass(),2) + 2 * (lepton.px() * pxNu + lepton.py() * pyNu);
275 double b1 = 2 * lepton.pz();
276 double A = 4*pow(lepton.E(),2) - b1*b1;
277 double B = -2 * c1 * b1;
278 double C = 4 * pow(lepton.E(), 2) * ptNu * ptNu - c1 * c1;
279 double discr = B*B - 4*A*C;
280 double r = 1;
281 double sol1, sol2;
282 if (discr > 0){
283 sol1 = (-B + sqrt(discr)) / (2*A);
284 sol2 = (-B - sqrt(discr)) / (2*A);
285 }
286 else {
287 // fitAlpha
288 std::valarray<double> pfix = { (m_mWpdg * m_mWpdg - lepton.mass() * lepton.mass()) / (2 * ptNu), met.phi() - lepton.phi(),
289 lepton.pt(), lepton.mass(), pxNu, pyNu };
290 MendelMin mm(delta2_fcn, 1, pfix);
291 mm.evolve(100);
292 valarray<double> fittest = mm.fittest();
293
294 const double alpha = fittest[0]*6.30-3.15; // map p[0] in [0,1] to alpha in [-3.15,3.15]
295 const double dphi = met.phi() - lepton.phi();
296 r = ( pow(m_mWpdg,2) - pow(lepton.mass(),2) );
297 r /= (2 * ptNu * (sqrt(pow(lepton.pt(),2) + pow(lepton.mass(),2)) - lepton.pt() * cos(dphi + alpha)));
298
299 const double old_p = ptNu;
300 const double old_phi = met.phi();
301 pxNu = r * old_p * cos(old_phi + alpha);
302 pyNu = r * old_p * sin(old_phi + alpha);
303 ptNu = sqrt (pxNu*pxNu + pyNu*pyNu);
304
305 c1 = pow(m_mWpdg,2) - pow(lepton.mass(),2) + 2 * (lepton.px() * pxNu + lepton.py() * pyNu);
306 B = -2 * c1 * b1;
307 C = 4 * pow(lepton.E(),2) * ptNu * ptNu - c1 * c1;
308 discr = B*B - 4*A*C;
309
310 sol1 = -B / (2*A);
311 sol2 = -B / (2*A);
312 }
313 // useSmallestPz
314 pzNu = (fabs(sol1) > fabs(sol2)) ? sol2 : sol1;
315
316 FourMomentum nu( sqrt(sqr(pxNu) + sqr(pyNu) + sqr(pzNu)), pxNu, pyNu, pzNu);
317
318 return nu;
319 }
320
321 };
322
323
324 RIVET_DECLARE_PLUGIN(ATLAS_2021_I1941095);
325
326}
|