1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6
  7#include "fastjet/JetDefinition.hh"
  8#include "fastjet/ClusterSequence.hh"
  9#include "fastjet/contrib/LundGenerator.hh"
 10
 11namespace Rivet {
 12
 13  /// @brief Lund jet plane with charged particles
 14  class ATLAS_2020_I1790256: public Analysis {
 15  public:
 16
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2020_I1790256);
 18
 19    /// @name Analysis methods
 20    //@{
 21
 22    void init() {
 23
 24      //Projections
 25      FinalState fs(Cuts::abseta < 4.5);
 26      FastJets jet4(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE);
 27      declare(jet4, "Jets");
 28
 29      ChargedFinalState tracks(Cuts::pT > 0.5*GeV && Cuts::abseta < 2.5);
 30      declare(tracks, "tracks");
 31
 32      book(_h_lundplane, 1,1,1);
 33
 34      _h_vs.resize(13);
 35      for (size_t i = 0; i < _h_vs.size(); ++i) {
 36        book(_h_vs[i] , i+3 , 1, 1);
 37      }
 38
 39      _h_hs.resize(19);
 40      for (size_t i = 0; i < _h_hs.size(); ++i) {
 41        book(_h_hs[i], i+16, 1, 1);
 42      }
 43
 44
 45      book(_njets, "_njets");
 46
 47    }
 48
 49    void analyze(const Event& event) {
 50
 51      const Jets jets = apply<JetAlg>(event, "Jets").jetsByPt(Cuts::pT > 300*GeV && Cuts::abseta < 2.1);
 52
 53      if (jets.size() < 2)  vetoEvent;
 54      if (jets[0].pT() < 675*GeV)  vetoEvent;
 55
 56      if ( (jets[0].pT()/jets[1].pT()) > 1.5 ) vetoEvent;
 57
 58       _njets->fill(2);
 59
 60       const Particles& tracks = apply<ChargedFinalState>(event, "tracks").particlesByPt();
 61
 62       Particles intracks1;
 63       Particles intracks2;
 64
 65       const Jet& j1 = jets[0];
 66       const Jet& j2 = jets[1];
 67
 68
 69      for (const Particle& p : tracks) {
 70        const double dr = deltaR(j1, p, PSEUDORAPIDITY);
 71        if (dr > 0.4) continue;
 72	      if (abs(p.pid()) == 13) continue;
 73        intracks1.push_back(p);
 74      }
 75
 76      for (const Particle& p : tracks) {
 77        const double dr = deltaR(j2, p, PSEUDORAPIDITY);
 78        if (dr > 0.4) continue;
 79        if (abs(p.pid()) == 13) continue;
 80        intracks2.push_back(p);
 81      }
 82
 83      JetDefinition tjet1_def(fastjet::cambridge_algorithm, 10);
 84      ClusterSequence tjet1_cs(intracks1, tjet1_def);
 85      vector<PseudoJet> tjets1 = fastjet::sorted_by_pt(tjet1_cs.inclusive_jets(0.0));
 86
 87      JetDefinition tjet2_def(fastjet::cambridge_algorithm, 10);
 88      ClusterSequence tjet2_cs(intracks2, tjet2_def);
 89      vector<PseudoJet> tjets2 = fastjet::sorted_by_pt(tjet2_cs.inclusive_jets(0.0));
 90
 91      if (tjets1.size() < 1 || tjets2.size() < 1) vetoEvent;
 92
 93      fjcontrib::LundGenerator lund;
 94      vector<fjcontrib::LundDeclustering> declusts1 = lund(tjets1[0]);
 95      for (size_t idecl = 0; idecl < declusts1.size(); ++idecl) {
 96        pair<double,double> coords = declusts1[idecl].lund_coordinates();
 97        double X = -0.9163 + coords.first;
 98        double Y = - log(declusts1[idecl].z());
 99
100        if (X > 0 && X < 4.33 && Y > log(1/0.5)  && Y < 8.6*log(1/0.5) ){
101
102          _h_lundplane->fill(X, Y);
103
104          double hdiv = (double)4.33/(double)13;
105          size_t i = floor(X/hdiv);
106          _h_vs[i]->fill(Y);
107
108          double vdiv = (8.6*log(1/0.5) - log(1/0.5))/(double)19;
109          size_t j = floor((Y - log(1/0.5))/vdiv);
110          _h_hs[j]->fill(X);
111
112        }
113      }
114
115      vector<fjcontrib::LundDeclustering> declusts2 = lund(tjets2[0]);
116      for (size_t idecl = 0; idecl < declusts2.size(); ++idecl) {
117        pair<double,double> coords = declusts2[idecl].lund_coordinates();
118
119        double X = -0.9163 + coords.first;
120        double Y = - log(declusts2[idecl].z());
121
122        if (X > 0 && X < 4.33 && Y > log(1/0.5)  && Y < 8.6*log(1/0.5) ) {
123
124          _h_lundplane->fill(X, Y);
125
126          double hdiv = (double)4.33/(double)13;
127          size_t i = floor(X/hdiv);
128          _h_vs[i]->fill(Y);
129
130          double vdiv = (8.6*log(1/0.5) - log(1/0.5))/(double)19;
131          size_t j = floor((Y - log(1/0.5))/vdiv);
132          _h_hs[j]->fill(X);
133        }
134      }
135    }
136
137
138    void finalize() {
139
140      double area = _njets->sumW();
141      if (area > 0) {
142        scale(_h_lundplane, 1/area);
143        scale(_h_vs, 1/(area*0.333));
144        scale(_h_hs, 1/(area*0.277));
145      }
146
147    }
148
149  private:
150
151
152    Histo2DPtr _h_lundplane;
153    vector<Histo1DPtr> _h_vs, _h_hs;
154    CounterPtr _njets;
155  };
156
157
158  // The hook for the plugin system
159  RIVET_DECLARE_PLUGIN(ATLAS_2020_I1790256);
160}