Rivet analyses referenceCMS_2021_I1901295Differential cross sections for top quark pair production in the full kinematic range using the lepton+jets final state in proton proton collisions at 13 TeVExperiment: CMS (LHC) Inspire ID: 1901295 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Measurements of differential and double-differential cross sections of top quark pair ($\mathrm{t\bar{t}}$) production are presented in the lepton+jets channels with a single electron or muon and jets in the final state. The analysis combines for the first time signatures of top quarks with low transverse momentum $p_\mathrm{T}$, where the top quark decay products can be identified as separated jets and isolated leptons, and with high $p_\mathrm{T}$, where the decay products are collimated and overlap. The measurements are based on proton-proton collision data at $\sqrt{s} = 13$ TeV collected by the CMS experiment at the LHC, corresponding to an integrated luminosity of 137 fb$^{-1}$. The cross sections are presented at the parton and particle levels, where the latter minimizes extrapolations based on theoretical assumptions. Most of the measured differential cross sections are well described by standard model predictions with the exception of some double-differential distributions. The inclusive $\mathrm{t\bar{t}}$ production cross section is measured to be $\sigma_\mathrm{t\bar{t}} = 791\pm25$ pb, which constitutes the most precise measurement in the lepton+jets channel to date. Source code: CMS_2021_I1901295.cc 1#include "Rivet/Analysis.hh"
2#include "Rivet/Math/LorentzTrans.hh"
3#include "Rivet/Projections/LeptonFinder.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/FinalState.hh"
6#include "Rivet/Projections/InvisibleFinalState.hh"
7#include "Rivet/Projections/PromptFinalState.hh"
8#include "Rivet/Tools/HistoGroup.hh"
9
10namespace Rivet {
11
12
13 /// ttbar lepton+jets at 13 TeV
14 class CMS_2021_I1901295 : public Analysis {
15 public:
16
17 RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2021_I1901295);
18
19 const double mtop_s = 172.5;
20 const double mw_s = 80.4;
21
22 void init() {
23
24 const FinalState fs(Cuts::abseta < 6.);
25
26 declare(fs, "FS");
27
28 declare(InvisibleFinalState(), "Invisibles");
29
30 FinalState all_photons(Cuts::abspid == PID::PHOTON);
31
32 PromptFinalState prompt_leptons(Cuts::abseta < 6. && (Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON));
33 prompt_leptons.acceptMuonDecays(true);
34 prompt_leptons.acceptTauDecays(true);
35
36 LeptonFinder dressed_leptons(prompt_leptons, all_photons, 0.1, Cuts::abseta < 2.4 && Cuts::pT > 15*GeV);
37 declare(dressed_leptons, "MyLeptons");
38
39 declare(FastJets(fs, JetAlg::ANTIKT, 0.4), "JetsAK4");
40 declare(FastJets(fs, JetAlg::ANTIKT, 0.8), "JetsAK8");
41
42 book(_h["thadpt"], "d159-x01-y01");
43 book(_h["tleppt"], "d163-x01-y01");
44 book(_h["thardpt"], "d167-x01-y01");
45 book(_h["tsoftpt"], "d171-x01-y01");
46 book(_h["st"], "d175-x01-y01");
47 book(_h["thady"], "d179-x01-y01");
48 book(_h["tlepy"], "d183-x01-y01");
49 book(_h["dy"], "d187-x01-y01");
50 book(_h["ady"], "d191-x01-y01");
51 book(_h["ttm"], "d195-x01-y01");
52 book(_h["ttpt"], "d199-x01-y01");
53 book(_h["tty"], "d203-x01-y01");
54 book(_h["ttphi"], "d207-x01-y01");
55 book(_h["cts"], "d211-x01-y01");
56 book(_h["lpt"], "d317-x01-y01");
57 book(_h["njet"], "d321-x01-y01");
58 book(_h["ht"], "d325-x01-y01");
59 book(_h["mevt"], "d329-x01-y01");
60
61 book(_h["thadpt_norm"], "d161-x01-y01");
62 book(_h["tleppt_norm"], "d165-x01-y01");
63 book(_h["thardpt_norm"], "d169-x01-y01");
64 book(_h["tsoftpt_norm"], "d173-x01-y01");
65 book(_h["st_norm"], "d177-x01-y01");
66 book(_h["thady_norm"], "d181-x01-y01");
67 book(_h["tlepy_norm"], "d185-x01-y01");
68 book(_h["dy_norm"], "d189-x01-y01");
69 book(_h["ady_norm"], "d193-x01-y01");
70 book(_h["ttm_norm"], "d197-x01-y01");
71 book(_h["ttpt_norm"], "d201-x01-y01");
72 book(_h["tty_norm"], "d205-x01-y01");
73 book(_h["ttphi_norm"], "d209-x01-y01");
74 book(_h["cts_norm"], "d213-x01-y01");
75 book(_h["lpt_norm"], "d319-x01-y01");
76 book(_h["njet_norm"], "d323-x01-y01");
77 book(_h["ht_norm"], "d327-x01-y01");
78 book(_h["mevt_norm"], "d331-x01-y01");
79
80 vector<double> thadptbins = {0.0, 80.0, 160.0, 240.0, 320.0, 400.0, 1500.0};
81 book(_b["thadpt_thady"], thadptbins);
82 book(_b["thadpt_thady_norm"], thadptbins);
83 for (size_t i = 0; i < _b["thadpt_thady"]->numBins(); ++i) {
84 book(_b["thadpt_thady"]->bin(i+1), 215 + i, 1, 1);
85 book(_b["thadpt_thady_norm"]->bin(i+1), 222 + i, 1, 1);
86 }
87 vector<double> ttmbins = {250.0, 420.0, 520.0, 620.0, 800.0, 1000.0, 3500.0};
88 book(_b["ttm_tty"], ttmbins);
89 book(_b["ttm_tty_norm"], ttmbins);
90 book(_b["ttm_cts"], ttmbins);
91 book(_b["ttm_cts_norm"], ttmbins);
92 book(_b["ttm_thadpt"], ttmbins);
93 book(_b["ttm_thadpt_norm"], ttmbins);
94 for (size_t i = 0; i < _b["ttm_tty"]->numBins(); ++i) {
95 book(_b["ttm_tty"]->bin(i+1), 229 + i, 1, 1);
96 book(_b["ttm_tty_norm"]->bin(i+1), 236 + i, 1, 1);
97 book(_b["ttm_cts"]->bin(i+1), 243 + i, 1, 1);
98 book(_b["ttm_cts_norm"]->bin(i+1), 250 + i, 1, 1);
99 book(_b["ttm_thadpt"]->bin(i+1), 257 + i, 1, 1);
100 book(_b["ttm_thadpt_norm"]->bin(i+1), 264 + i, 1, 1);
101 }
102 vector<double> ttptbins = {0.0, 50.0, 120.0, 200.0, 300.0, 400.0, 1200.0};
103 book(_b["ttpt_thadpt"], ttptbins);
104 book(_b["ttpt_thadpt_norm"], ttptbins);
105 for (size_t i = 0; i < _b["ttpt_thadpt"]->numBins(); ++i) {
106 book(_b["ttpt_thadpt"]->bin(i+1), 271 + i, 1, 1);
107 book(_b["ttpt_thadpt_norm"]->bin(i+1), 278 + i, 1, 1);
108 }
109 vector<double> adybins = {0.0, 0.6, 1.2, 1.8, 3.5};
110 book(_b["ady_ttm"], adybins);
111 book(_b["ady_ttm_norm"], adybins);
112 for (size_t i = 0; i < _b["ady_ttm"]->numBins(); ++i) {
113 book(_b["ady_ttm"]->bin(i+1), 285 + i, 1, 1);
114 book(_b["ady_ttm_norm"]->bin(i+1), 290 + i, 1, 1);
115 }
116 vector<double> ttmbins2 = {250.0, 700.0, 1000.0, 3500.0};
117 book(_b["ttm_dy"], ttmbins2);
118 book(_b["ttm_dy_norm"], ttmbins2);
119 for (size_t i = 0; i < _b["ttm_dy"]->numBins(); ++i) {
120 book(_b["ttm_dy"]->bin(i+1), 295 + i, 1, 1);
121 book(_b["ttm_dy_norm"]->bin(i+1), 299 + i, 1, 1);
122 }
123 vector<double> topybins = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5};
124 book(_b["topy_topbary"], topybins);
125 book(_b["topy_topbary_norm"], topybins);
126 for (size_t i = 0; i < _b["topy_topbary"]->numBins(); ++i) {
127 book(_b["topy_topbary"]->bin(i+1), 303 + i, 1, 1);
128 book(_b["topy_topbary_norm"]->bin(i+1), 311 + i, 1, 1);
129 }
130 vector<double> njetbins = {-0.5, 0.5, 1.5, 2.5, 3.5};
131 book(_b["njet_ttm"], njetbins);
132 book(_b["njet_ttm_norm"], njetbins);
133 book(_b["njet_ttpt"], njetbins);
134 book(_b["njet_ttpt_norm"], njetbins);
135 book(_b["njet_thadpt"], njetbins);
136 book(_b["njet_thadpt_norm"], njetbins);
137 for (size_t i = 0; i < _b["njet_ttm"]->numBins(); ++i) {
138 book(_b["njet_ttm"]->bin(i+1), 333 + i, 1, 1);
139 book(_b["njet_ttm_norm"]->bin(i+1), 338 + i, 1, 1);
140 book(_b["njet_ttpt"]->bin(i+1), 343 + i, 1, 1);
141 book(_b["njet_ttpt_norm"]->bin(i+1), 348 + i, 1, 1);
142 book(_b["njet_thadpt"]->bin(i+1), 353 + i, 1, 1);
143 book(_b["njet_thadpt_norm"]->bin(i+1), 358 + i, 1, 1);
144 }
145 }
146
147 void analyze(const Event &event) {
148 DressedLeptons m_leptons;
149 Jets m_bjets;
150 Jets m_ljets;
151 Jets m_alljets;
152 Particles m_thadboosted;
153 Particles m_tlepboosted;
154 Particles m_additionaljets;
155
156 int numvetoleps = 0;
157 const DressedLeptons &dressedleptons =
158 apply<LeptonFinder>(event, "MyLeptons").dressedLeptons();
159 for (const DressedLepton &lep : dressedleptons) {
160 if (lep.pt()/GeV > 30. && lep.abseta() < 2.4) {
161 m_leptons.push_back(lep);
162 } else {
163 ++numvetoleps;
164 }
165 }
166 if (m_leptons.size() != 1 || numvetoleps != 0) {
167 vetoEvent;
168 }
169
170 DressedLepton lepton = m_leptons[0];
171
172 const Particles &invfspars =
173 apply<FinalState>(event, "Invisibles").particles(Cuts::abseta < 6.);
174 FourMomentum nusum = accumulate(
175 invfspars.begin(), invfspars.end(), FourMomentum(0., 0., 0., 0.),
176 [&](const FourMomentum &invmom, const Particle &par) {
177 return invmom + par.momentum();
178 });
179
180 const Jets &allAK4Jets =
181 apply<FastJets>(event, "JetsAK4")
182 .jetsByPt(Cuts::abseta < 2.4 && Cuts::pT > 25. * GeV);
183
184 for (const Jet &jet : allAK4Jets) {
185 if (deltaR(lepton, jet) > 0.4) {
186 if (jet.bTagged()) {
187 m_bjets.push_back(jet);
188 } else {
189 m_ljets.push_back(jet);
190 }
191
192 m_alljets.push_back(jet);
193 } else if (jet.bTagged() && lepton.pt()/GeV > 50.) {
194 const Particle &undressedlep(lepton.bareLepton());
195 bool injet = false;
196 for (const Particle &con : jet.particles()) {
197 if (con.momentum() == undressedlep.momentum()) {
198 injet = true;
199 break;
200 }
201 }
202 FourMomentum purejet = jet.momentum();
203 if (injet) {
204 purejet -= lepton.momentum();
205 }
206 double checkreco;
207 FourMomentum nureco =
208 numom(nusum, lepton.momentum(), purejet, checkreco);
209 if (checkreco < 0.) {
210 continue;
211 }
212 FourMomentum tlepmom(nureco + purejet + lepton.momentum());
213 if (!injet) {
214 tlepmom += lepton.momentum();
215 }
216 if (tlepmom.pt()/GeV > 380. && tlepmom.mass()/GeV > 100. &&
217 tlepmom.mass()/GeV < 240.) {
218 m_tlepboosted.emplace_back(6, tlepmom);
219 }
220 }
221 }
222
223 const Jets &allAK8Jets =
224 apply<FastJets>(event, "JetsAK8")
225 .jetsByPt(Cuts::abseta < 2.4 && Cuts::pT > 380. * GeV);
226 for (const Jet &jet : allAK8Jets) {
227 if (deltaR(lepton, jet) < 0.8) {
228 continue;
229 }
230 if (jet.bTagged() && jet.mass()/GeV > 120.) {
231 m_thadboosted.emplace_back(6, jet.momentum());
232 }
233 }
234
235 bool reco = false;
236 Particle thad;
237 Particle tlep;
238
239 // boosted thad, boosted tlep reconstruction
240 if (m_tlepboosted.size() != 0) {
241 double bestdm = numeric_limits<double>::max();
242 sort(m_tlepboosted.begin(), m_tlepboosted.end(),
243 [&mtop = mtop_s](const FourMomentum &A, const FourMomentum &B) {
244 return abs(A.mass()/GeV - mtop) < abs(B.mass()/GeV - mtop);
245 });
246 tlep = m_tlepboosted[0];
247 for (const Particle &th : m_thadboosted) {
248 if (deltaR(th, tlep) < 1.2) {
249 continue;
250 }
251 double dm = abs(th.mass()/GeV - mtop_s);
252 if (dm < bestdm) {
253 bestdm = dm;
254 thad = Particle(6, th.momentum());
255 reco = true;
256 }
257 }
258 }
259
260 if (reco) {
261 for (const Jet &jet : m_alljets) {
262 if (deltaR(thad, jet) < 1.2) {
263 continue;
264 }
265 if (deltaR(tlep, jet) < 0.6) {
266 continue;
267 }
268 m_additionaljets.emplace_back(1, jet.momentum());
269 }
270 } else {
271 // boosted thad, resolved tlep reconstruction
272 double bestcoma = numeric_limits<double>::max();
273 Particle tlepcoma;
274 Particle thadcoma;
275 Particle blepcoma;
276 for (const Particle &th : m_thadboosted) {
277 for (const Jet &bjl : m_bjets) {
278 if (deltaR(th, bjl) < 1.2) {
279 continue;
280 }
281 if (deltaR(th, lepton) < 1.2) {
282 continue;
283 }
284 double checkreco;
285 FourMomentum nureco = numom(nusum, lepton.momentum(), bjl, checkreco);
286 if (checkreco < 0.) {
287 continue;
288 }
289 Particle tl =
290 Particle(6, nureco + lepton.momentum() + bjl.momentum());
291 if (tl.mass()/GeV < 100. || tl.mass()/GeV > 240.) {
292 continue;
293 }
294 double coma = pow(th.mass()/GeV - mtop_s, 2) + pow(tl.mass()/GeV - mtop_s, 2);
295 if (coma < bestcoma) {
296 bestcoma = coma;
297 blepcoma = Particle(5, bjl.momentum());
298 //!!!To reproduce the distributions in the paper, we have to use
299 //!nusum as neutrino momentum here instead of nureco.!!!
300 tlepcoma = Particle(6, nusum + lepton.momentum() + bjl.momentum());
301 // tlepcoma = Particle(6, tl);
302 thadcoma = th;
303 }
304 }
305 }
306
307 double bestcomb = numeric_limits<double>::max();
308 Particle tlepcomb;
309 Particle thadcomb;
310 Particles ttdecay(4);
311 if (m_bjets.size() >= 2 && m_ljets.size() >= 2) {
312 // resolved thad, resolved tlep reconstruction
313 for (const Jet &bjl : m_bjets) {
314 double checkreco;
315 FourMomentum nureco =
316 numom(nusum, lepton.momentum(), bjl.momentum(), checkreco);
317 if (checkreco < 0.) {
318 continue;
319 }
320 FourMomentum tl(lepton.momentum() + nureco + bjl.momentum());
321 if (tl.mass()/GeV < 100. || tl.mass()/GeV > 240.) {
322 continue;
323 }
324
325 for (size_t a = 0; a < m_ljets.size(); ++a) {
326 const Jet &lja = m_ljets[a];
327 for (size_t b = 0; b < a; ++b) {
328 const Jet &ljb = m_ljets[b];
329 FourMomentum wh(lja.momentum() + ljb.momentum());
330 for (const Jet &bjh : m_bjets) {
331 if (&bjh == &bjl) {
332 continue;
333 }
334 FourMomentum th(wh + bjh.momentum());
335 if (th.mass()/GeV < 100. || th.mass()/GeV > 240.) {
336 continue;
337 }
338
339 double comb = pow(wh.mass()/GeV - mw_s, 2) +
340 pow(th.mass()/GeV - mtop_s, 2) +
341 pow(tl.mass()/GeV - mtop_s, 2);
342 if (comb < bestcomb) {
343 bestcomb = comb;
344 thadcomb = Particle(6, th);
345 //!!!To reproduce the distributions in the paper, we have to
346 //!use nusum as neutrino momentum here instead of nureco.!!!
347 tlepcomb =
348 Particle(6, lepton.momentum() + bjl.momentum() + nusum);
349 // tlepcomb = Particle(6, tl);
350 ttdecay[0] = Particle(5, bjh);
351 ttdecay[1] = Particle(1, lja);
352 ttdecay[2] = Particle(1, ljb);
353 ttdecay[3] = Particle(5, bjl);
354 }
355 }
356 }
357 }
358 }
359 }
360
361 if (bestcoma != numeric_limits<double>::max() &&
362 bestcomb != numeric_limits<double>::max()) {
363 if (abs(thadcoma.mass()/GeV - mtop_s) < abs(thadcomb.mass()/GeV - mtop_s)) {
364 bestcomb = numeric_limits<double>::max();
365 } else {
366 bestcoma = numeric_limits<double>::max();
367 }
368 }
369
370 if (bestcoma != numeric_limits<double>::max()) {
371 reco = true;
372 thad = thadcoma;
373 tlep = tlepcoma;
374 for (const Jet &jet : m_alljets) {
375 if (deltaR(blepcoma, jet) < 0.01) {
376 continue;
377 }
378 if (deltaR(thadcoma, jet) < 1.2) {
379 continue;
380 }
381 m_additionaljets.emplace_back(1, jet.momentum());
382 }
383 } else if (bestcomb != numeric_limits<double>::max()) {
384 reco = true;
385 thad = thadcomb;
386 tlep = tlepcomb;
387 for (const Jet &jet : m_alljets) {
388 if (find_if(ttdecay.begin(), ttdecay.end(), [&](const Particle &par) {
389 return deltaR(jet, par) < 0.01;
390 }) != ttdecay.end()) {
391 continue;
392 }
393 m_additionaljets.emplace_back(1, jet.momentum());
394 }
395 }
396 }
397
398 if (!reco) {
399 vetoEvent;
400 }
401 FourMomentum tt(thad.momentum() + tlep.momentum());
402 FourMomentum top = lepton.pid() > 0 ? thad : tlep;
403 FourMomentum topbar = lepton.pid() < 0 ? thad : tlep;
404
405 const double ht = sum(m_additionaljets, Kin::pT, 0.)/GeV;
406 FourMomentum mevt = std::accumulate(m_additionaljets.begin(), m_additionaljets.end(), tt,
407 [](const FourMomentum &mevt, const Particle &jet) {
408 return mevt + jet;
409 });
410
411 LorentzTransform boostcms;
412 boostcms.setBetaVec(-1. * tt.betaVec());
413 FourMomentum topcms = boostcms(top);
414 double cts = topcms.vector3().dot(tt.vector3()) / (topcms.p() * tt.p());
415 double day = abs(top.rapidity()) - abs(topbar.rapidity());
416 double ady = abs(top.rapidity() - topbar.rapidity());
417 double ttphi = abs(deltaPhi(top, topbar)) * 180. / M_PI;
418
419 dualfill("thadpt", thad.pt()/GeV);
420 dualfill("tleppt", tlep.pt()/GeV);
421 dualfill("thardpt", max(thad.pt()/GeV, tlep.pt()/GeV));
422 dualfill("tsoftpt", min(thad.pt()/GeV, tlep.pt()/GeV));
423 dualfill("st", thad.pt()/GeV + tlep.pt()/GeV);
424 dualfill("thady", abs(thad.rapidity()));
425 dualfill("tlepy", abs(tlep.rapidity()));
426 dualfill("dy", day);
427 dualfill("ady", ady);
428 dualfill("ttm", tt.mass()/GeV);
429 dualfill("ttpt", tt.pt()/GeV);
430 dualfill("tty", abs(tt.rapidity()));
431 dualfill("ttphi", ttphi);
432 dualfill("cts", cts);
433 dualfill("lpt", lepton.pt()/GeV);
434 dualfill("njet", min(m_additionaljets.size(), 6u) + 0.5);
435 dualfill("ht", ht/GeV);
436 dualfill("mevt", mevt.mass()/GeV);
437
438 dualfill("thadpt_thady", thad.pt()/GeV, abs(thad.rapidity()));
439 dualfill("ttm_tty", tt.mass()/GeV, abs(tt.rapidity()));
440 dualfill("ttm_cts", tt.mass()/GeV, cts);
441 dualfill("ttm_thadpt", tt.mass()/GeV, thad.pt()/GeV);
442 dualfill("ttpt_thadpt", tt.pt()/GeV, thad.pt()/GeV);
443 dualfill("ady_ttm", ady, tt.mass()/GeV);
444 dualfill("ttm_dy", tt.mass()/GeV, day);
445 dualfill("topy_topbary", abs(top.rapidity()), abs(topbar.rapidity()));
446 dualfill("njet_ttm", min(m_additionaljets.size(), 3u), tt.mass()/GeV);
447 dualfill("njet_ttpt", min(m_additionaljets.size(), 3u), tt.pt()/GeV);
448 dualfill("njet_thadpt", min(m_additionaljets.size(), 3u), thad.pt()/GeV);
449 }
450
451 void finalize() {
452
453 const double sf = crossSection()/picobarn/sumOfWeights();
454
455 for (auto& item : _h) {
456 if (item.first.find("_norm") != string::npos) {
457 normalize(item.second, 1.0, false);
458 }
459 else scale(item.second, sf);
460 }
461
462 for (auto& item : _b) {
463 if (item.first.find("_norm") != string::npos) {
464 normalizeGroup(item.second, 1.0, false);
465 }
466 else {
467 scale(item.second, sf);
468 }
469 divByGroupWidth(item.second);
470 }
471 }
472
473 map<string,Histo1DPtr> _h;
474 map<string,Histo1DGroupPtr> _b;
475
476 void dualfill(const string& tag, const double value) {
477 _h[tag]->fill(value); _h[tag + "_norm"]->fill(value);
478 }
479
480 void dualfill(const string& tag, const double val1, const double val2) {
481 _b[tag]->fill(val1, val2); _b[tag + "_norm"]->fill(val1, val2);
482 }
483
484 FourMomentum numom(const FourMomentum &met, const FourMomentum &l,
485 const FourMomentum &b, double &mtres) {
486 mtres = -1.;
487 FourMomentum va;
488 FourMomentum vb;
489
490 double tltn = l.px() * met.px() + l.py() * met.py();
491 double tntn = met.px() * met.px() + met.py() * met.py();
492
493 double C = pow(l.pz() / l.E(), 2) - 1.;
494 double B =
495 pow(mw_s / l.E(), 2) * l.pz() + 2. * l.pz() / l.E() / l.E() * tltn;
496 double A = -1. * tntn + pow(mw_s * mw_s / 2. / l.E(), 2) +
497 pow(mw_s / l.E(), 2) * tltn + pow(tltn / l.E(), 2);
498
499 double D = pow(B / C / 2., 2) - A / C;
500 if (D >= 0.) {
501 va.setXYZM(met.px(), met.py(), -0.5 * B / C + sqrt(D), 0.);
502 vb.setXYZM(met.px(), met.py(), -0.5 * B / C - sqrt(D), 0.);
503 double ma = (va + l + b).mass()/GeV;
504 double mb = (vb + l + b).mass()/GeV;
505
506 if (abs(ma - mtop_s) < abs(mb - mtop_s)) {
507 mtres = ma;
508 return va;
509 } else {
510 mtres = mb;
511 return vb;
512 }
513 } else {
514 double As = 0.25 * (pow(mw_s * mw_s * l.pz() / l.E() / l.E(), 2) / C -
515 pow(mw_s * mw_s / l.E(), 2));
516 double Bsx =
517 (pow(mw_s * l.pz() / l.E() / l.E(), 2) / C - pow(mw_s / l.E(), 2)) *
518 l.px() * met.px();
519 double Bsy =
520 (pow(mw_s * l.pz() / l.E() / l.E(), 2) / C - pow(mw_s / l.E(), 2)) *
521 l.py() * met.py();
522 double Csxx = (pow(l.pz() / l.E() / l.E(), 2) / C - 1. / l.E() / l.E()) *
523 l.px() * met.px() * l.px() * met.px() +
524 met.px() * met.px();
525 double Csxy = (pow(l.pz() / l.E() / l.E(), 2) / C - 1. / l.E() / l.E()) *
526 2. * l.px() * met.px() * l.py() * met.py();
527 double Csyy = (pow(l.pz() / l.E() / l.E(), 2) / C - 1. / l.E() / l.E()) *
528 l.py() * met.py() * l.py() * met.py() +
529 met.py() * met.py();
530
531 As /= C;
532 Bsx /= C;
533 Bsy /= C;
534 Csxx /= C;
535 Csxy /= C;
536 Csyy /= C;
537
538 double x = 1.;
539 double y = 1.;
540
541 double U =
542 As + x * Bsx + y * Bsy + x * x * Csxx + x * y * Csxy + y * y * Csyy;
543
544 double step = 0.1;
545
546 while (true) {
547 double dx = Bsx + 2. * Csxx * x + Csxy * y;
548 double dy = Bsy + 2. * Csyy * y + Csxy * x;
549
550 x += step * dx / sqrt(dx * dx + dy * dy);
551 y += step * dy / sqrt(dx * dx + dy * dy);
552
553 double nU = (As + x * Bsx + y * Bsy + x * x * Csxx + x * y * Csxy +
554 y * y * Csyy);
555
556 if (nU * U < 0.) {
557 step *= -0.5;
558 }
559 U = nU;
560 if (abs(U) < 0.01) {
561 break;
562 }
563 }
564 double pz = -0.5 / C *
565 (pow(mw_s / l.E(), 2) * l.pz() +
566 2. * l.pz() / l.E() / l.E() *
567 (l.px() * met.px() * x + l.py() * met.py() * y));
568 va.setXYZM(met.px() * x, met.py() * y, pz, 0.);
569 mtres = (va + l + b).mass()/GeV;
570 return va;
571 }
572
573 return va;
574 }
575
576 };
577
578 RIVET_DECLARE_PLUGIN(CMS_2021_I1901295);
579
580} // namespace Rivet
|