1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// Jet fragmentation at 7 TeV
 11  class ATLAS_2011_I929691 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I929691);
 16
 17
 18    /// Initialisation
 19    void init() {
 20      const FinalState fs(Cuts::abseta < 2.0);
 21
 22      FastJets antikt_06_jets(fs, JetAlg::ANTIKT, 0.6, JetMuons::NONE, JetInvisibles::NONE);
 23      declare(antikt_06_jets, "jets");
 24
 25      ChargedFinalState tracks(Cuts::pT > 0.5*GeV && Cuts::abseta < 2.0);
 26      declare(tracks, "tracks");
 27
 28      // Set up the histograms (each element is a binning in jet pT)
 29      for (size_t i = 0; i < 10; ++i) {
 30        book(_p_F_z[i]    , i+ 1, 1, 1);
 31        book(_p_rho_r[i]  , i+11, 1, 1);
 32        book(_p_f_pTrel[i], i+21, 1, 1);
 33      }
 34
 35    }
 36
 37
 38    // Per-event analysis
 39    void analyze(const Event& event) {
 40
 41      const Jets alljets = apply<FastJets>(event, "jets").jetsByPt(Cuts::absrap < 1.2);
 42      const Particles& tracks = apply<ChargedFinalState>(event, "tracks").particlesByPt();
 43
 44      for (size_t i = 0; i < 10; ++i) {
 45
 46        const Jets jets = select(alljets, Cuts::pT > bedges[i] && Cuts::pT < bedges[i+1]);
 47        const int n_jets = jets.size();
 48        if (n_jets == 0) continue;
 49
 50        // First... count the tracks
 51        Histo1D h_ntracks_z(_p_F_z[i]->binning()),
 52                h_ntracks_r(_p_rho_r[i]->binning()),
 53                h_ntracks_pTrel(_p_f_pTrel[i]->binning());
 54        for (const Jet& j : jets) {
 55          for (const Particle& p : tracks) {
 56            const double dr = deltaR(j, p, RAPIDITY);
 57            if (dr > 0.6) continue; // The paper uses pseudorapidity, but this is a requirement for filling the histogram
 58            h_ntracks_z.fill(z(j, p), 1.0/n_jets);
 59            h_ntracks_r.fill(dr, 1.0/n_jets);
 60            h_ntracks_pTrel.fill(pTrel(j, p), 1.0/n_jets);
 61          }
 62        }
 63
 64        // Then... calculate the observable and fill the profiles
 65        for (const auto& b : h_ntracks_z.bins())
 66          _p_F_z[i]->fill(b.xMid(), b.sumW()/b.dVol());
 67        for (const auto& b : h_ntracks_r.bins())
 68          _p_rho_r[i]->fill(b.xMid(), b.sumW()/annulus_area(b.xMin(), b.xMax()));
 69        for (const auto& b : h_ntracks_pTrel.bins())
 70          _p_f_pTrel[i]->fill(b.xMid(), b.sumW()/b.dVol());
 71
 72      }
 73
 74    }
 75
 76
 77    double z (const Jet& jet, const Particle& ch) {
 78      return dot(jet.p3(), ch.p3()) / jet.p3().mod2();
 79    }
 80
 81    double pTrel (const Jet& jet, const Particle& ch) {
 82      return (ch.p3().cross(jet.p3())).mod()/(jet.p3().mod());
 83    }
 84
 85    // To calculate the area of the annulus in an r bin
 86    double annulus_area(double r1, double r2) {
 87      return M_PI*(sqr(r2) - sqr(r1));
 88    }
 89
 90
 91  private:
 92
 93    Profile1DPtr _p_F_z[10], _p_rho_r[10], _p_f_pTrel[10];
 94    const vector<double> bedges = { 25., 40., 60., 80., 110., 160., 210., 260., 310., 400., 500. };
 95
 96  };
 97
 98
 99  RIVET_DECLARE_PLUGIN(ATLAS_2011_I929691);
100
101
102}