Rivet analyses referenceATLAS_2019_I1724098Jet substructure at 13 TeVExperiment: ATLAS (LHC) Inspire ID: 1724098 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
A measurement of jet substructure observables is presented using data collected in 2016 by the ATLAS experiment at the LHC with proton-proton collisions at $\sqrt{s} = 13$ TeV. Large-radius jets groomed with the trimming and soft-drop algorithms are studied. Dedicated event selections are used to study jets produced by light quarks or gluons, and hadronically decaying top quarks and W bosons. The observables measured are sensitive to substructure, and therefore are typically used for tagging large-radius jets from boosted massive particles. These include the energy correlation functions and the $N$-subjettiness variables. The number of subjets and the Les Houches angularity are also considered. The distributions of the substructure variables, corrected for detector effects, are compared to the predictions of various Monte Carlo event generators. They are also compared between the large-radius jets originating from light quarks or gluons, and hadronically decaying top quarks and $W$ bosons. Source code: ATLAS_2019_I1724098.cc 1// -*- C++ -*-
2
3#include "Rivet/Analysis.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/PromptFinalState.hh"
6#include "Rivet/Projections/ChargedLeptons.hh"
7#include "Rivet/Projections/DressedLeptons.hh"
8#include "Rivet/Projections/MissingMomentum.hh"
9
10#include "fastjet/tools/Filter.hh"
11#include "fastjet/contrib/SoftDrop.hh"
12#include "fastjet/contrib/Nsubjettiness.hh"
13#include "fastjet/contrib/Njettiness.hh"
14#include "fastjet/contrib/EnergyCorrelator.hh"
15
16
17namespace Rivet {
18
19
20 /// @brief Jet substructure at 13 TeV
21 class ATLAS_2019_I1724098: public Analysis {
22 public:
23
24 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1724098);
25
26 private:
27
28 /// @name Analysis methods
29 //@{
30
31 void init() {
32
33 _mode = 0; // default is do eveything
34 if ( getOption("MODE") == "DJ" ) _mode = 1;
35 else if ( getOption("MODE") == "TW" ) _mode = 2;
36
37 //Projections
38 const FinalState fs(Cuts::abseta < 4.5);
39
40 // Get photons used to dress leptons
41 const FinalState photons(Cuts::abspid == PID::PHOTON);
42
43 // Use all bare muons as input to the DressedMuons projection
44 PromptFinalState bare_mu(Cuts::abspid == PID::MUON, true);
45 PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, true);
46
47 // Muons must have |eta| < 2.5
48 Cut eta_ranges = Cuts::abseta < 2.5;
49 DressedLeptons dressed_mu(photons, bare_mu, 0.1, eta_ranges && Cuts::pT > 30*GeV, true);
50 declare(dressed_mu, "muons");
51 DressedLeptons dressed_el(photons, bare_el, 0.1, eta_ranges && Cuts::pT > 25*GeV, true);
52 declare(dressed_el, "electrons");
53
54 FastJets fj(fs, FastJets::ANTIKT, 1.0, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE);
55 declare(fj, "FJets");
56 FastJets sj(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE);
57 declare(sj, "Jets");
58
59 ChargedLeptons lfs(FinalState(Cuts::abseta < 2.5 && Cuts::pT > 25*GeV));
60 declare(lfs, "LFS");
61
62 MissingMomentum missmom(fs);
63 declare(missmom, "MissingMomentum");
64
65 _trimmer = fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.05));
66
67 // Dijet
68 if (_mode == 0 || _mode == 1) {
69 // SD
70 book(_h["dj_sdnsj"],1,1,1);
71 book(_h["dj_sdlha"] ,2,1,1);
72 book(_h["dj_sdc2"] ,3,1,1);
73 book(_h["dj_sdd2"] ,4,1,1);
74 book(_h["dj_sdecf2"] ,5,1,1);
75 book(_h["dj_sdecf3"] ,6,1,1);
76 // TRIMMED
77 book(_h["dj_nsj"] ,23,1,1);
78 book(_h["dj_lha"] ,24,1,1);
79 book(_h["dj_c2"] ,25,1,1);
80 book(_h["dj_d2"] ,26,1,1);
81 book(_h["dj_ecf2"] ,27,1,1);
82 book(_h["dj_ecf3"] ,28,1,1);
83 }
84
85 if (_mode == 0 || _mode == 2) {
86 //Top SD
87 book(_h["tw_sdnsj"] ,7,1,1);
88 book(_h["tw_sdlha"] ,8,1,1);
89 book(_h["tw_sdc2"] ,9,1,1);
90 book(_h["tw_sdd2"] ,10,1,1);
91 book(_h["tw_sdecf2"] ,11,1,1);
92 book(_h["tw_sdecf3"] ,12,1,1);
93 book(_h["tw_sdtau21"] ,13,1,1);
94 book(_h["tw_sdtau32"] ,14,1,1);
95 //W SD
96 book(_h["tw_wsdnsj"] ,15,1,1);
97 book(_h["tw_wsdlha"] ,16,1,1);
98 book(_h["tw_wsdc2"] ,17,1,1);
99 book(_h["tw_wsdd2"] ,18,1,1);
100 book(_h["tw_wsdecf2"] ,19,1,1);
101 book(_h["tw_wsdecf3"] ,20,1,1);
102 book(_h["tw_wsdtau21"] ,21,1,1);
103 book(_h["tw_wsdtau32"] ,22,1,1);
104 //Top TRIMMED
105 book(_h["tw_nsj"] ,29,1,1);
106 book(_h["tw_lha"] ,30,1,1);
107 book(_h["tw_c2"] ,31,1,1);
108 book(_h["tw_d2"] ,32,1,1);
109 book(_h["tw_ecf2"] ,33,1,1);
110 book(_h["tw_ecf3"] ,34,1,1);
111 book(_h["tw_tau21"] ,35,1,1);
112 book(_h["tw_tau32"] ,36,1,1);
113 //W TRIMMED
114 book(_h["tw_wnsj"] ,37,1,1);
115 book(_h["tw_wlha"] ,38,1,1);
116 book(_h["tw_wc2"] ,39,1,1);
117 book(_h["tw_wd2"] ,40,1,1);
118 book(_h["tw_wecf2"] ,41,1,1);
119 book(_h["tw_wecf3"] ,42,1,1);
120 book(_h["tw_wtau21"] ,43,1,1);
121 book(_h["tw_wtau32"] ,44,1,1);
122 }
123
124 }
125
126 /// Do the analysis
127 void analyze(const Event& event) {
128
129 if (_mode == 0 || _mode == 1) doDIJET(event);
130 if (_mode == 0 || _mode == 2) doTW(event);
131
132 }
133
134
135 void doDIJET(const Event& event) {
136
137 double nsub, lha, ecf2, ecf3, c2, d2;
138 double sdnsub, sdlha, sdecf2, sdecf3, sdc2, sdd2;
139
140 lha = sdlha = 0.0;
141
142 const double beta = 1;
143
144 const Particles & leptons = apply<ChargedLeptons>(event, "LFS").particles();
145 if (leptons.size()) return;
146
147 // Normal fatjets
148 const Jets &fjets = apply<JetAlg>(event, "FJets").jetsByPt();
149
150 // Trim the fatjets
151 PseudoJets tr_ljets;
152 for (size_t i = 0; i < fjets.size(); ++i) {
153 tr_ljets += _trimmer(fjets[i]);
154 tr_ljets[tr_ljets.size()-1].set_user_index(i);
155 }
156
157 size_t nBaseline = count(tr_ljets, [](const Jet &j) { return j.pT() > 200*GeV && j.abseta() < 2.5; });
158 if (nBaseline < 2) return;
159
160 ifilter_select(tr_ljets, [](const PseudoJet &j) { return j.perp() > 450*GeV; });
161 if (tr_ljets.size() > 1) tr_ljets = sorted_by_pt(tr_ljets);
162 else if (tr_ljets.empty()) return;
163
164 if (abs(tr_ljets[0].eta()) > 1.5) return;
165 const fastjet::PseudoJet &LJet = tr_ljets[0];
166 size_t uindex = tr_ljets[0].user_index();
167
168 // Nsubjets
169 JetDefinition subjet_def(fastjet::kt_algorithm, 0.2);
170 ClusterSequence subjet_cs(LJet.constituents(), subjet_def);
171 PseudoJets subjets = sorted_by_pt(subjet_cs.inclusive_jets(10.0));
172 nsub = subjets.size();
173
174 // LHA
175 for (const PseudoJet& p : LJet.constituents()) {
176 double trpt = p.pt();
177 double trtheta = p.squared_distance(LJet);
178 lha += pow(trpt, 1.0) * pow(trtheta, 0.25);
179 }
180 double lterm = pow(LJet.pt(), 1.0) * pow(1.0, 0.5);
181 if (lterm !=0) lha /= lterm;
182 else lha = -99;
183
184 // C2
185 fastjet::contrib::EnergyCorrelator ECF3(3,beta,fastjet::contrib::EnergyCorrelator::pt_R);
186 fastjet::contrib::EnergyCorrelator ECF2(2,beta,fastjet::contrib::EnergyCorrelator::pt_R);
187 fastjet::contrib::EnergyCorrelator ECF1(1,beta,fastjet::contrib::EnergyCorrelator::pt_R);
188 double recf3 = ECF3(LJet);
189 double recf2 = ECF2(LJet);
190 double recf1 = ECF1(LJet);
191 c2 = (recf2 != 0 ? recf3 * recf1 / (recf2*recf2) : -1);
192 d2 = (recf2 != 0 ? recf3 * (recf1*recf1*recf1) /(recf2*recf2*recf2) : -1);
193 ecf2 = (recf1 != 0 ? recf2 / (recf1*recf1) : -1);
194 ecf3 = (recf1 != 0 ? recf3 / (recf1*recf1*recf1) : -1);
195
196 // Fill Histograms for trimmed
197 _h["dj_nsj"]->fill(nsub);
198 _h["dj_c2"]->fill(c2);
199 _h["dj_d2"]->fill(d2);
200 _h["dj_lha"]->fill(lha);
201 _h["dj_ecf2"]->fill(ecf2);
202 _h["dj_ecf3"]->fill(ecf3);
203
204
205 ////////////////////////////////////////////
206 // Soft Drop
207 fastjet::contrib::SoftDrop sd(0.0, 0.1);
208 PseudoJet SDLJet = sd(fjets[uindex]);
209
210 ClusterSequence subjet_sdcs(SDLJet.constituents(), subjet_def);
211
212 PseudoJets sdsubjets = sorted_by_pt(subjet_sdcs.inclusive_jets(10.0));
213 sdnsub = sdsubjets.size();
214
215 for (const PseudoJet& sd_p : SDLJet.constituents()) {
216 double spt = sd_p.pt();
217 double stheta = sd_p.squared_distance(SDLJet);
218 sdlha += pow(spt, 1.0) * pow(stheta, 0.25);
219 }
220
221 double sdterm = pow(SDLJet.pt(), 1.0) * pow(1.0, 0.5);
222 if (sdterm !=0) sdlha /= sdterm;
223 else sdlha = -99;
224
225 double sdrecf3 = ECF3(SDLJet);
226 double sdrecf2 = ECF2(SDLJet);
227 double sdrecf1 = ECF1(SDLJet);
228
229 sdc2 = (sdrecf2 != 0 ? sdrecf3 * sdrecf1 / (sdrecf2*sdrecf2) : -1);
230 sdd2 = (sdrecf2 != 0 ? sdrecf3 * (sdrecf1*sdrecf1*sdrecf1) /(sdrecf2*sdrecf2*sdrecf2) : -1);
231
232 sdecf2 = (sdrecf1 !=0 ? sdrecf2 /(sdrecf1*sdrecf1) : -1);
233 sdecf3 = (sdrecf1 !=0 ? sdrecf3 / (sdrecf1*sdrecf1*sdrecf1) : -1);
234
235 _h["dj_sdnsj"]->fill(sdnsub);
236 _h["dj_sdc2"]->fill(sdc2);
237 _h["dj_sdd2"]->fill(sdd2);
238 _h["dj_sdlha"]->fill(sdlha);
239 _h["dj_sdecf2"]->fill(sdecf2);
240 _h["dj_sdecf3"]->fill(sdecf3);
241
242 }
243
244 void doTW(const Event& event) {
245
246 double nsub, lha, tau21, tau32, ecf2, ecf3, c2, d2;
247 double sdnsub, sdlha, sdtau21, sdtau32, sdecf2, sdecf3, sdc2, sdd2;
248 double wnsub, wlha, wtau21, wtau32, wecf2, wecf3, wc2, wd2;
249 double wsdnsub, wsdlha, wsdtau21, wsdtau32, wsdecf2, wsdecf3, wsdc2, wsdd2;
250
251 lha = sdlha = wlha = wsdlha = 0.0;
252
253 const double beta = 1, Rcut = 1;
254
255 const vector<DressedLepton>& muons = apply<DressedLeptons>(event, "muons").dressedLeptons();
256 if (muons.size() != 1) return;
257
258 const vector<DressedLepton>& electrons = apply<DressedLeptons>(event, "electrons").dressedLeptons();
259 if (electrons.size() != 0) return;
260
261 const FourMomentum muonmom = muons[0].momentum();
262 const MissingMomentum& missmom = apply<MissingMomentum>(event, "MissingMomentum");
263 FourMomentum missvec = missmom.missingMomentum();
264 double met = missmom.missingPt();
265 if (met < 20*GeV) return;
266 const double transmass = sqrt( 2 * muons[0].pT() * met * (1 - cos(deltaPhi(muons[0], missvec))) );
267 if (transmass + met <= 60*GeV) return;
268
269 const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
270 if(jets.empty()) return;
271
272 int lepJetIndex = -1;
273 for (size_t i = 0; i < jets.size(); ++i) {
274 const Jet& jet = jets[i];
275 if ((deltaR(jet, muons[0]) < 1.5) && (deltaR(jet, muons[0]) > 0.4) ) {
276 lepJetIndex = i;
277 break;
278 }
279 }
280 if (lepJetIndex < 0) return;
281 const Jet& lepjet = jets[lepJetIndex];
282
283
284 const Jets fjets = apply<JetAlg>(event, "FJets").jetsByPt();
285 PseudoJets tr_ljets_all;
286 for (const Jet& j : fjets) {
287 tr_ljets_all += _trimmer(j);
288 }
289
290 PseudoJets tr_ljets;
291 for (size_t i = 0; i < tr_ljets_all.size(); ++i) {
292 const PseudoJet tj = tr_ljets_all[i];
293 if (tj.perp() > 150*GeV && fabs(tj.eta()) < 2.5) {
294 tr_ljets += tj;
295 tr_ljets[tr_ljets.size()-1].set_user_index(i);
296 }
297 }
298
299 if (tr_ljets.size() < 1) return;
300 if (tr_ljets.size() > 1) tr_ljets = sorted_by_pt(tr_ljets);
301
302 const Jet& fjet = tr_ljets[0];
303 size_t uindex = tr_ljets[0].user_index();
304
305 const double dR_fatjet = deltaR(lepjet, fjet);
306 const double dPhi_fatjet = deltaPhi(muons[0], fjet);
307
308 if (dR_fatjet < 1.5 || dPhi_fatjet < 2.3) return;
309
310 Jets bjets, non_bjets;
311 for (const Jet& jet : jets)
312 (jet.bTagged() ? bjets : non_bjets) += jet;
313 if (bjets.empty()) return;
314
315 double min_bdR = 99;
316 int bindex = 0;
317 int k=0;
318
319 for (const Jet& bjet : bjets) {
320 double bdR = deltaR(fjet, bjet);
321 if(bdR < min_bdR){
322 min_bdR = bdR;
323 bindex = k;
324 }
325 k++;
326 }
327
328 size_t tw = 0;
329 double dR_fjet_bjet = deltaR(fjet, bjets[bindex]);
330 if (fjet.abseta() < 1.5) {
331 if (dR_fjet_bjet < 1.0 && fjet.mass() > 140 && fjet.pT() > 350*GeV) tw = 1;
332 if (dR_fjet_bjet > 1.0 && dR_fjet_bjet < 1.8 && fjet.mass() < 100*GeV && fjet.mass() > 60*GeV && fjet.pT() > 200*GeV) tw = 2;
333 }
334
335
336 // Top plots:
337 if (tw==1) {
338
339 PseudoJet LJet = fjet;
340 JetDefinition subjet_def(fastjet::kt_algorithm, 0.2);
341 ClusterSequence subjet_cs(LJet.constituents(), subjet_def);
342 PseudoJets subjets = sorted_by_pt(subjet_cs.inclusive_jets(10.0));
343 nsub = subjets.size();
344
345 // LHA
346 for (const PseudoJet& p : LJet.constituents()){
347 double pt = p.pt();
348 double theta = p.squared_distance(LJet);
349 lha += pow(pt, 1.0) * pow(theta, 0.25);
350 }
351 double lterm = pow(LJet.pt(), 1.0) * pow(1.0, 0.5);
352 if (lterm) lha /= lterm;
353 else lha = -99;
354
355
356 // NSubjettiness
357 fastjet::contrib::Nsubjettiness nSub1(1, fastjet::contrib::OnePass_WTA_KT_Axes(), fastjet::contrib::NormalizedMeasure(beta,Rcut));
358 fastjet::contrib::Nsubjettiness nSub2(2, fastjet::contrib::OnePass_WTA_KT_Axes(), fastjet::contrib::NormalizedMeasure(beta,Rcut));
359 fastjet::contrib::Nsubjettiness nSub3(3, fastjet::contrib::OnePass_WTA_KT_Axes(), fastjet::contrib::NormalizedMeasure(beta,Rcut));
360 double tau1 = nSub1.result(LJet);
361 double tau2 = nSub2.result(LJet);
362 double tau3 = nSub3.result(LJet);
363 if(tau1 != 0) tau21 = tau2/tau1;
364 else tau21 = -99;
365 if(tau2 != 0) tau32 = tau3/tau2;
366 else tau32 = -99;
367
368
369 //C2
370 fastjet::contrib::EnergyCorrelator ECF3(3,beta,fastjet::contrib::EnergyCorrelator::pt_R);
371 fastjet::contrib::EnergyCorrelator ECF2(2,beta,fastjet::contrib::EnergyCorrelator::pt_R);
372 fastjet::contrib::EnergyCorrelator ECF1(1,beta,fastjet::contrib::EnergyCorrelator::pt_R);
373
374 double recf3 = ECF3(LJet);
375 double recf2 = ECF2(LJet);
376 double recf1 = ECF1(LJet);
377
378
379 c2 = (recf2 != 0 ? recf3 * recf1 / (recf2*recf2) : -1);
380 d2 = (recf2 != 0 ? recf3 * (recf1*recf1*recf1) /(recf2*recf2*recf2) : -1);
381
382 ecf2 = (recf1 !=0 ? recf2 /(recf1*recf1) : -1);
383 ecf3 = (recf1 !=0 ? recf3 / (recf1*recf1*recf1) : -1);
384
385 _h["tw_nsj"]->fill(nsub);
386 _h["tw_lha"]->fill(lha);
387 _h["tw_tau21"]->fill(tau21);
388 _h["tw_tau32"]->fill(tau32);
389 _h["tw_c2"]->fill(c2);
390 _h["tw_d2"]->fill(d2);
391 _h["tw_ecf2"]->fill(ecf2);
392 _h["tw_ecf3"]->fill(ecf3);
393
394
395 // Soft Drop
396 fastjet::contrib::SoftDrop sd(0.0, 0.1);
397 PseudoJet SDLJet = sd(fjets[uindex]);
398 ClusterSequence subjet_sdcs(SDLJet.constituents(), subjet_def);
399
400 PseudoJets sdsubjets = sorted_by_pt(subjet_sdcs.inclusive_jets(10.0));
401 sdnsub = sdsubjets.size();
402
403 for (const PseudoJet& sd_p : SDLJet.constituents()){
404 double spt = sd_p.pt();
405 double stheta = sd_p.squared_distance(SDLJet);
406 sdlha += pow(spt, 1.0) * pow(stheta, 0.25);
407 }
408
409 double sdlterm = pow(SDLJet.pt(), 1.0) * pow(1.0, 0.5);
410 if (sdlterm) sdlha /= sdlterm;
411 else sdlha = -99;
412
413
414 double sdtau1 = nSub1.result(SDLJet);
415 double sdtau2 = nSub2.result(SDLJet);
416 double sdtau3 = nSub3.result(SDLJet);
417 if(sdtau1 != 0) sdtau21 = sdtau2/sdtau1;
418 else sdtau21 = -99;
419 if(sdtau2 != 0) sdtau32 = sdtau3/sdtau2;
420 else sdtau32 = -99;
421
422 double sdrecf3 = ECF3(SDLJet);
423 double sdrecf2 = ECF2(SDLJet);
424 double sdrecf1 = ECF1(SDLJet);
425
426 sdc2 = (sdrecf2 != 0 ? sdrecf3 * sdrecf1 / (sdrecf2*sdrecf2) : -1);
427 sdd2 = (sdrecf2 != 0 ? sdrecf3 * (sdrecf1*sdrecf1*sdrecf1) /(sdrecf2*sdrecf2*sdrecf2) : -1);
428
429 sdecf2 = (sdrecf1 !=0 ? sdrecf2 /(sdrecf1*sdrecf1) : -1);
430 sdecf3 = (sdrecf1 !=0 ? sdrecf3 / (sdrecf1*sdrecf1*sdrecf1) : -1);
431
432 _h["tw_sdnsj"]->fill(sdnsub);
433 _h["tw_sdlha"]->fill(sdlha);
434 _h["tw_sdtau21"]->fill(sdtau21);
435 _h["tw_sdtau32"]->fill(sdtau32);
436 _h["tw_sdc2"]->fill(sdc2);
437 _h["tw_sdd2"]->fill(sdd2);
438 _h["tw_sdecf2"]->fill(sdecf2);
439 _h["tw_sdecf3"]->fill(sdecf3);
440
441 }
442
443
444 // W plots
445
446 if(tw ==2){
447
448 PseudoJet LJet = fjet;
449 JetDefinition subjet_def(fastjet::kt_algorithm, 0.2);
450 ClusterSequence subjet_cs(LJet.constituents(), subjet_def);
451
452 PseudoJets subjets = sorted_by_pt(subjet_cs.inclusive_jets(10.0));
453 wnsub = subjets.size();
454
455 // LHA
456 for (const PseudoJet& wp : LJet.constituents()){
457 double wpt = wp.pt();
458 double wtheta = wp.squared_distance(fjet);
459 wlha += pow(wpt, 1.0) * pow(wtheta, 0.25);
460 }
461
462 double fterm = pow(fjet.pt(), 1.0) * pow(1.0, 0.5);
463 if (fterm) wlha /= fterm;
464 else wlha = -99;
465
466
467
468 // NSubjettiness
469
470 fastjet::contrib::Nsubjettiness nSub1(1, fastjet::contrib::OnePass_WTA_KT_Axes(), fastjet::contrib::NormalizedMeasure(beta,Rcut));
471 fastjet::contrib::Nsubjettiness nSub2(2, fastjet::contrib::OnePass_WTA_KT_Axes(), fastjet::contrib::NormalizedMeasure(beta,Rcut));
472 fastjet::contrib::Nsubjettiness nSub3(3, fastjet::contrib::OnePass_WTA_KT_Axes(), fastjet::contrib::NormalizedMeasure(beta,Rcut));
473 double wtau1 = nSub1.result(LJet);
474 double wtau2 = nSub2.result(LJet);
475 double wtau3 = nSub3.result(LJet);
476 if(wtau1 != 0) wtau21 = wtau2/wtau1;
477 else wtau21 = -99;
478 if(wtau2 != 0) wtau32 = wtau3/wtau2;
479 else wtau32 = -99;
480
481 //C2
482 fastjet::contrib::EnergyCorrelator ECF3(3,beta,fastjet::contrib::EnergyCorrelator::pt_R);
483 fastjet::contrib::EnergyCorrelator ECF2(2,beta,fastjet::contrib::EnergyCorrelator::pt_R);
484 fastjet::contrib::EnergyCorrelator ECF1(1,beta,fastjet::contrib::EnergyCorrelator::pt_R);
485
486 double wrecf3 = ECF3(LJet);
487 double wrecf2 = ECF2(LJet);
488 double wrecf1 = ECF1(LJet);
489
490 wc2 = (wrecf2 != 0 ? wrecf3 * wrecf1 / (wrecf2*wrecf2) : -1);
491 wd2 = (wrecf2 != 0 ? wrecf3 * (wrecf1*wrecf1*wrecf1) /(wrecf2*wrecf2*wrecf2) : -1);
492
493 wecf2 = (wrecf1 !=0 ? wrecf2 /(wrecf1*wrecf1) : -1);
494 wecf3 = (wrecf1 !=0 ? wrecf3 / (wrecf1*wrecf1*wrecf1) : -1);
495
496 _h["tw_wnsj"]->fill(wnsub);
497 _h["tw_wlha"]->fill(wlha);
498 _h["tw_wtau21"]->fill(wtau21);
499 _h["tw_wtau32"]->fill(wtau32);
500 _h["tw_wc2"]->fill(wc2);
501 _h["tw_wd2"]->fill(wd2);
502 _h["tw_wecf2"]->fill(wecf2);
503 _h["tw_wecf3"]->fill(wecf3);
504
505
506 //SD
507 fastjet::contrib::SoftDrop sd(0.0, 0.1);
508 PseudoJet SDLJet = sd(fjets[uindex]);
509 ClusterSequence subjet_sdcs(SDLJet.constituents(), subjet_def);
510
511 PseudoJets sdsubjets = sorted_by_pt(subjet_sdcs.inclusive_jets(10.0));
512 wsdnsub = sdsubjets.size();
513
514 for (const PseudoJet& sd_p : SDLJet.constituents()){
515 double spt = sd_p.pt();
516 double stheta = sd_p.squared_distance(SDLJet);
517 wsdlha += pow(spt, 1.0) * pow(stheta, 0.25);
518 }
519
520 double wsdlterm = pow(SDLJet.pt(), 1.0) * pow(1.0, 0.5);
521 if (wsdlterm) wsdlha /= wsdlterm;
522 else wsdlha = -99;
523
524
525 double wsdtau1 = nSub1.result(SDLJet);
526 double wsdtau2 = nSub2.result(SDLJet);
527 double wsdtau3 = nSub3.result(SDLJet);
528 if (wsdtau1 != 0) wsdtau21 = wsdtau2/wsdtau1;
529 else wsdtau21 = -99;
530 if (wsdtau2 != 0) wsdtau32 = wsdtau3/wsdtau2;
531 else wsdtau32 = -99;
532
533 double wsdrecf3 = ECF3(SDLJet);
534 double wsdrecf2 = ECF2(SDLJet);
535 double wsdrecf1 = ECF1(SDLJet);
536
537 wsdc2 = (wsdrecf2 != 0 ? wsdrecf3 * wsdrecf1 / (wsdrecf2*wsdrecf2) : -1);
538 wsdd2 = (wsdrecf2 != 0 ? wsdrecf3 * (wsdrecf1*wsdrecf1*wsdrecf1) /(wsdrecf2*wsdrecf2*wsdrecf2) : -1);
539
540 wsdecf2 = (wsdrecf1 !=0 ? wsdrecf2 /(wsdrecf1*wsdrecf1) : -1);
541 wsdecf3 = (wsdrecf1 !=0 ? wsdrecf3 / (wsdrecf1*wsdrecf1*wsdrecf1) : -1);
542
543 _h["tw_wsdnsj"]->fill(wsdnsub);
544 _h["tw_wsdlha"]->fill(wsdlha);
545 _h["tw_wsdtau21"]->fill(wsdtau21);
546 _h["tw_wsdtau32"]->fill(wsdtau32);
547 _h["tw_wsdc2"]->fill(wsdc2);
548 _h["tw_wsdd2"]->fill(wsdd2);
549 _h["tw_wsdecf2"]->fill(wsdecf2);
550 _h["tw_wsdecf3"]->fill(wsdecf3);
551
552 }
553 }
554
555
556 void finalize() {
557 for (auto hist : _h) { normalize(hist.second); }
558 }
559
560 private:
561
562 fastjet::Filter _trimmer;
563 map<string, Histo1DPtr> _h;
564
565 protected:
566 size_t _mode;
567 };
568
569 // The hook for the plugin system
570 RIVET_DECLARE_PLUGIN(ATLAS_2019_I1724098);
571}
|