1// -*- C++ -*-
  2
  3#include "Rivet/Analysis.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief ATLAS inclusive jet pT spectrum, di-jet Mass and di-jet chi
 10  class ATLAS_2010_I871366: public Analysis {
 11  public:
 12
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2010_I871366);
 14
 15    /// Choice of jet algorithm
 16    enum Alg { AKT4=0, AKT6=1 };
 17
 18
 19    void init() {
 20      FinalState fs;
 21      declare(fs, "FinalState");
 22      declare(FastJets(fs, JetAlg::ANTIKT, 0.6), "AntiKT06");
 23      declare(FastJets(fs, JetAlg::ANTIKT, 0.4), "AntiKT04");
 24
 25      const vector<double> ybins{ 0.0, 0.3, 0.8, 1.2, 2.1, 2.8 };
 26      const vector<double> massBinsForChi{ 340, 520, 800, 1200 };
 27
 28
 29      size_t ptDsOffset(0), massDsOffset(10), chiDsOffset(20);
 30      for (size_t alg = 0; alg < 2; ++alg) {
 31        book(_pThistos[alg], ybins);
 32        for (auto& b : _pThistos[alg]->bins()) {
 33          book(b, b.index() + ptDsOffset, 1, 1);
 34        }
 35        ptDsOffset += 5;
 36
 37        book(_massVsY[alg], ybins);
 38        for (auto& b : _massVsY[alg]->bins()) {
 39          book(b, b.index() + massDsOffset, 1, 1);
 40        }
 41        massDsOffset += 5;
 42
 43        book(_chiVsMass[alg], massBinsForChi);
 44        for (auto& b : _chiVsMass[alg]->bins()) {
 45          book(b, b.index() + chiDsOffset, 1, 1);
 46        }
 47        chiDsOffset += 3;
 48      }
 49    }
 50
 51
 52    void analyze(const Event& evt) {
 53      Jets jetAr[2];
 54      jetAr[AKT6] = apply<FastJets>(evt, "AntiKT06").jetsByPt(Cuts::pT > 30*GeV);
 55      jetAr[AKT4] = apply<FastJets>(evt, "AntiKT04").jetsByPt(Cuts::pT > 30*GeV);
 56
 57      // Identify the dijets
 58      for (size_t alg = 0; alg < 2; ++alg) {
 59        vector<FourMomentum> leadjets;
 60        for (const Jet& jet : jetAr[alg]) {
 61          const double pT = jet.pT();
 62          const double absy = jet.absrap();
 63          _pThistos[alg]->fill(absy, pT/GeV);
 64
 65          if (absy < 2.8 && leadjets.size() < 2) {
 66            if (leadjets.empty() && pT < 60*GeV) continue;
 67            leadjets.push_back(jet.momentum());
 68          }
 69        }
 70
 71        // Veto if no acceptable dijet found
 72        if (leadjets.size() < 2) {
 73          MSG_DEBUG("Could not find two suitable leading jets");
 74          continue;
 75        }
 76
 77        const double rap1 = leadjets[0].rapidity();
 78        const double rap2 = leadjets[1].rapidity();
 79        const double mass = (leadjets[0] + leadjets[1]).mass();
 80        const double ymax = max(fabs(rap1), fabs(rap2));
 81        const double chi = exp(fabs(rap1 - rap2));
 82        if (fabs(rap1 + rap2) < 2.2) {
 83          _chiVsMass[alg]->fill(mass/GeV, chi);
 84        }
 85        _massVsY[alg]->fill(ymax, mass/GeV);
 86
 87      }
 88    }
 89
 90
 91    void finalize() {
 92      const double sf = crossSectionPerEvent()/picobarn;
 93      // factor 0.5 needed because it is differential in dy and not d|y|
 94      scale(_pThistos, 0.5*sf);
 95      scale(_massVsY, sf);
 96      scale(_chiVsMass, sf);
 97
 98      divByGroupWidth(_pThistos);
 99      divByGroupWidth(_massVsY);
100      divByGroupWidth(_chiVsMass);
101    }
102
103
104  private:
105
106    /// The inclusive pT spectrum for akt6 and akt4 jets (array index is jet type from enum above)
107    Histo1DGroupPtr _pThistos[2];
108
109    /// The di-jet mass spectrum binned in rapidity for akt6 and akt4 jets (array index is jet type from enum above)
110    Histo1DGroupPtr _massVsY[2];
111
112    /// The di-jet chi distribution binned in mass for akt6 and akt4 jets (array index is jet type from enum above)
113    Histo1DGroupPtr _chiVsMass[2];
114
115  };
116
117
118
119  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2010_I871366, ATLAS_2010_S8817804);
120
121}