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