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