1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  class ATLAS_2012_I1082936 : public Analysis {
 10  public:
 11
 12    /// @name Constructors etc.
 13    /// @{
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1082936);
 17
 18    /// @}
 19
 20
 21  public:
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    /// Book histograms and initialise projections before the run
 27    void init() {
 28
 29      const FinalState fs;
 30      declare(fs,"FinalState");
 31
 32      FastJets fj04(fs,  JetAlg::ANTIKT, 0.4);
 33      fj04.useInvisibles();
 34      declare(fj04, "AntiKT04");
 35
 36      FastJets fj06(fs,  JetAlg::ANTIKT, 0.6);
 37      fj06.useInvisibles();
 38      declare(fj06, "AntiKT06");
 39
 40
 41      // Histogram booking copied from the previous analysis
 42      const vector<double> ybins{ 0.0, 0.3, 0.8, 1.2, 2.1, 2.8, 3.6, 4.4 };
 43      const vector<double> ystarbins{ 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.4};
 44
 45      size_t ptDsOffset(0), massDsOffset(2);
 46      for (size_t alg = 0; alg < 2; ++alg) {
 47        book(_pThistos[alg], ybins);
 48        for (auto& b : _pThistos[alg]->bins()) {
 49          book(b, 1+ptDsOffset, 1, b.index());
 50        }
 51        ++ptDsOffset;
 52
 53        book(_mass[alg], ystarbins);
 54        for (auto& b : _mass[alg]->bins()) {
 55          book(b, 1+massDsOffset, 1, b.index());
 56        }
 57        ++massDsOffset;
 58      }
 59    }
 60
 61    /// Perform the per-event analysis
 62    void analyze(const Event& event) {
 63
 64      Jets jetAr[2];
 65      jetAr[AKT6] = apply<FastJets>(event, "AntiKT06").jetsByPt(Cuts::pT > 20*GeV);
 66      jetAr[AKT4] = apply<FastJets>(event, "AntiKT04").jetsByPt(Cuts::pT > 20*GeV);
 67
 68      // Loop over jet "radii" used in analysis
 69      for (size_t alg = 0; alg < 2; ++alg) {
 70        // Identify dijets
 71        vector<FourMomentum> leadjets;
 72        for (const Jet& jet : jetAr[alg]) {
 73          const double pT = jet.pT();
 74          const double absy = jet.absrap();
 75          _pThistos[alg]->fill(absy, pT/GeV);
 76
 77          if (absy < 4.4 && leadjets.size() < 2) {
 78            if (leadjets.empty() && pT < 30*GeV) continue;
 79            leadjets.push_back(jet.momentum());
 80          }
 81        }
 82        // Make sure we have a leading jet with pT >30 GeV and a second to leading jet with pT>20 GeV
 83        if (leadjets.size() < 2) {
 84          MSG_DEBUG("Could not find two suitable leading jets");
 85          continue;
 86        }
 87
 88        const double y1 = leadjets[0].rapidity();
 89        const double y2 = leadjets[1].rapidity();
 90        const double ystar = fabs(y1-y2)/2.;
 91        const double m = (leadjets[0] + leadjets[1]).mass();
 92        // Fill mass histogram
 93        _mass[alg]->fill(ystar, m/TeV);
 94      }
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      // factor 0.5 needed because it is differential in dy and not d|y|
101      scale(_pThistos, 0.5*crossSectionPerEvent()/picobarn);
102      scale(_mass, crossSectionPerEvent()/picobarn);
103      divByGroupWidth(_pThistos);
104      divByGroupWidth(_mass);
105    }
106
107    /// @}
108
109
110  private:
111
112    // Data members like post-cuts event weight counters go here
113
114    enum Alg { AKT4=0, AKT6=1 };
115
116  private:
117
118    /// The inclusive pT spectrum for akt6 and akt4 jets (array index is jet type from enum above)
119    Histo1DGroupPtr _pThistos[2];
120
121    /// The di-jet mass spectrum binned in rapidity for akt6 and akt4 jets (array index is jet type from enum above)
122    Histo1DGroupPtr _mass[2];
123  };
124
125  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1082936);
126
127}