rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2019_I1724098

Jet substructure at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1724098
Status: VALIDATED
Authors:
  • Deepak Kar
  • Amal Vaidya
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • Multijet / ttbar production at 13 TeV

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}