1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/VisibleFinalState.hh"
  6#include "Rivet/Projections/JetShape.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief ATLAS jet shape analysis
 12  /// @author Andy Buckley, Judith Katzy, Francesc Vives
 13  class ATLAS_2011_I882984 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I882984);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    void init() {
 24      // Set up projections
 25      const FinalState fs((Cuts::etaIn(-5.0, 5.0)));
 26      declare(fs, "FS");
 27      FastJets fj(fs, JetAlg::ANTIKT, 0.6);
 28      fj.useInvisibles();
 29      declare(fj, "Jets");
 30
 31      // Specify pT bins
 32      _ptedges = {{ 30.0, 40.0, 60.0, 80.0, 110.0, 160.0, 210.0, 260.0, 310.0, 400.0, 500.0, 600.0 }};
 33      _yedges  = {{ 0.0, 0.3, 0.8, 1.2, 2.1, 2.8 }};
 34
 35      // Register a jet shape projection and histogram for each pT bin
 36      for (size_t i = 0; i < 11; ++i) {
 37        for (size_t j = 0; j < 6; ++j) {
 38          if (i == 8  && j == 4) continue;
 39          if (i == 9  && j == 4) continue;
 40          if (i == 10 && j != 5) continue;
 41
 42          // Set up projections for each (pT,y) bin
 43          _jsnames_pT[i][j] = "JetShape" + to_str(i) + to_str(j);
 44          const double ylow = (j < 5) ? _yedges[j] : _yedges.front();
 45          const double yhigh = (j < 5) ? _yedges[j+1] : _yedges.back();
 46          const JetShape jsp(fj, 0.0, 0.7, 7, _ptedges[i], _ptedges[i+1], ylow, yhigh, RAPIDITY);
 47          declare(jsp, _jsnames_pT[i][j]);
 48
 49          // Book profile histograms for each (pT,y) bin
 50          book(_profhistRho_pT[i][j] ,i+1, j+1, 1);
 51          book(_profhistPsi_pT[i][j] ,i+1, j+1, 2);
 52        }
 53      }
 54    }
 55
 56
 57
 58    /// Do the analysis
 59    void analyze(const Event& evt) {
 60
 61      // Get jets and require at least one to pass pT and y cuts
 62      const Jets jets = apply<FastJets>(evt, "Jets")
 63        .jetsByPt(Cuts::ptIn(_ptedges.front()*GeV, _ptedges.back()*GeV) && Cuts::absrap < 2.8);
 64      MSG_DEBUG("Jet multiplicity before cuts = " << jets.size());
 65      if (jets.size() == 0) {
 66        MSG_DEBUG("No jets found in required pT and rapidity range");
 67        vetoEvent;
 68      }
 69      // Calculate and histogram jet shapes
 70      for (size_t ipt = 0; ipt < 11; ++ipt) {
 71        for (size_t jy = 0; jy < 6; ++jy) {
 72          if (ipt == 8 && jy == 4) continue;
 73          if (ipt == 9 && jy == 4) continue;
 74          if (ipt == 10 && jy != 5) continue;
 75          const JetShape& jsipt = apply<JetShape>(evt, _jsnames_pT[ipt][jy]);
 76          for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) {
 77            for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
 78              const double r_rho = jsipt.rBinMid(rbin);
 79              _profhistRho_pT[ipt][jy]->fill(r_rho, (1./0.1)*jsipt.diffJetShape(ijet, rbin));
 80              const double r_Psi = jsipt.rBinMid(rbin);
 81              _profhistPsi_pT[ipt][jy]->fill(r_Psi, jsipt.intJetShape(ijet, rbin));
 82            }
 83          }
 84        }
 85      }
 86    }
 87
 88    /// @}
 89
 90
 91  private:
 92
 93    /// @name Analysis data
 94    /// @{
 95
 96    /// Jet \f$ p_\perp\f$ bins.
 97    vector<double> _ptedges; // This can't be a raw array if we want to initialise it non-painfully
 98    vector<double> _yedges;
 99
100    /// JetShape projection name for each \f$p_\perp\f$ bin.
101    string _jsnames_pT[11][6];
102
103    ///@}
104
105
106    /// @name Histograms
107    /// @{
108    Profile1DPtr _profhistRho_pT[11][6];
109    Profile1DPtr _profhistPsi_pT[11][6];
110    /// @}
111
112  };
113
114
115
116  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I882984, ATLAS_2011_S8924791);
117
118}