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