1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  namespace {
 10
 11    inline double calcz(const Jet& j, const Particle& p) {
 12      const double num = j.p3().dot(p.p3());
 13      const double den = j.p3().mod2();
 14      return num/den;
 15    }
 16
 17    inline double calcptrel(const Jet& j, const Particle& p) {
 18      const double num = j.p3().cross(p.p3()).mod();
 19      const double den = j.p3().mod();
 20      return num/den;
 21    }
 22
 23    inline double calcr(const Jet& j, const Particle& p) {
 24      return deltaR(j.rapidity(), j.phi(), p.rapidity(), p.phi());
 25    }
 26
 27    // For annulus area kludge
 28    /// @todo Improve somehow... need normalisation *without* bin width factors!
 29    inline double calcrweight(const Jet& j, const Particle& p) {
 30      size_t nBins_r = 26;
 31      double bins_r[] = { 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10,
 32                          0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24, 0.26, 0.28, 0.30,
 33                          0.35, 0.40, 0.45, 0.50, 0.55, 0.60 };
 34      double r = calcr(j,p);
 35      for (size_t bin = 0 ; bin < nBins_r ; bin++) {
 36        if (r < bins_r[bin+1]) {
 37          double up = bins_r[bin+1];
 38          double down = bins_r[bin];
 39          return ((up-down)/(M_PI*(up*up-down*down)));
 40        }
 41      }
 42      return 1.0;
 43    }
 44  }
 45
 46
 47  class ATLAS_2011_I919017 : public Analysis {
 48  public:
 49
 50    /// @name Constructors etc.
 51    /// @{
 52
 53    /// Constructor
 54    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I919017);
 55
 56    /// @}
 57
 58
 59  public:
 60
 61    /// @name Analysis methods
 62    /// @{
 63
 64    /// Book histograms and initialise projections before the run
 65    void init() {
 66      ChargedFinalState cfs((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  0.3*GeV));
 67      FastJets trkjets04(cfs, JetAlg::ANTIKT, 0.4);
 68      FastJets trkjets06(cfs, JetAlg::ANTIKT, 0.6);
 69      declare(trkjets04, "Jets04");
 70      declare(trkjets06, "Jets06");
 71
 72      // Book histograms
 73      const vector<string> rapbins{"00_05_", "05_10_", "10_15_", "15_19_"};
 74      const vector<string> ptbins{"04_06", "06_10", "10_15", "15_24", "24_40"};
 75      for (const string R : {"04_", "06_"}) {
 76        size_t iy = (R == "04_")? 1 : 2;
 77        book(_c[R+"sumw"], "/TMP/"+R+"sumw");
 78        size_t irap = 0;
 79        for (const string& rapbin : rapbins) {
 80          ++irap;
 81          book(_hinc[R+rapbin], irap, 1, iy);
 82          size_t ipt = 0;
 83          for (const string& ptbin : ptbins) {
 84            ++ipt;
 85            const string suff = R + rapbin + ptbin;
 86            if (irap == 1) {
 87              size_t offset = rapbins.size();
 88              book(_h["N"+R+"00_19_"+ptbin], offset+ipt, 1, iy);
 89              offset += ptbins.size()*ptbins.size();
 90              book(_h["z"+R+"00_19_"+ptbin], offset+ipt, 1, iy);
 91              offset += ptbins.size()*ptbins.size();
 92              book(_h["ptrel"+R+"00_19_"+ptbin], offset+ipt, 1, iy);
 93              offset += ptbins.size()*ptbins.size();
 94              book(_h["rdA"+R+"00_19_"+ptbin], offset+ipt, 1, iy);
 95              book(_c["numjets"+R+"00_19_"+ptbin], "/TMP/numjets"+R+"00_19_"+ptbin);
 96            }
 97            size_t offset = rapbins.size() + irap*ptbins.size();
 98            book(_h["N"+suff], offset+ipt, 1, iy);
 99            offset = rapbins.size() + (irap+ptbins.size())*ptbins.size();
100            book(_h["z"+suff], offset+ipt, 1, iy);
101            offset = rapbins.size() + (irap+2*ptbins.size())*ptbins.size();
102            book(_h["ptrel"+suff], offset+ipt, 1, iy);
103            offset = rapbins.size() + (irap+3*ptbins.size())*ptbins.size();
104            book(_h["rdA"+suff], offset+ipt, 1, iy);
105            book(_c["numjets"+suff], "/TMP/numjets"+suff);
106          }
107        }
108      }
109    }
110
111
112    /// Perform the per-event analysis
113    void analyze(const Event& event) {
114      doCollection("04", event);
115      doCollection("06", event);
116    }
117
118    void doCollection(const string& Rval, const Event& event) {
119      const Jets& jets = apply<JetFinder>(event, "Jets"+Rval).jets();
120      if (jets.empty())  return;
121      const string R = Rval+"_";
122      _c[R+"sumw"]->fill();
123      for (const Jet& j : jets) {
124        const double jetpt = j.pT()/GeV;
125        const double jetrap = j.absrap();
126        string rapbin, ptbin;
127        if (jetrap >= 1.9)  continue;
128
129        if (jetrap < 0.5)       rapbin = "00_05_";
130        else if (jetrap < 1.0)  rapbin = "05_10_";
131        else if (jetrap < 1.5)  rapbin = "10_15_";
132        else                    rapbin = "15_19_";
133        _hinc[R+rapbin]->fill(jetpt);
134
135        if (inRange(jetpt, 4., 6.))         ptbin = "04_06";
136        else if (inRange(jetpt, 6., 10.))   ptbin = "06_10";
137        else if (inRange(jetpt, 10., 15.))  ptbin = "10_15";
138        else if (inRange(jetpt, 15., 24.))  ptbin = "15_24";
139        else if (inRange(jetpt, 24., 40.))  ptbin = "24_40";
140        else  continue;
141        fillgroup(j, R+"00_19_"+ptbin);
142        fillgroup(j, R+rapbin+ptbin);
143      } // each jet
144    }
145
146    void fillgroup(const Jet& j, const string& suff) {
147      _c["numjets"+suff]->fill();
148      _h["N"+suff]->fill(j.particles().size());
149      for (const Particle& p : j.particles()) {
150        _h["z"+suff]->fill(calcz(j, p));
151        _h["ptrel"+suff]->fill(calcptrel(j, p));
152        _h["rdA"+suff]->fill(calcr(j, p), calcrweight(j, p));
153      }
154    }
155
156    /// Normalise histograms etc., after the run
157    void finalize() {
158
159      const double xsec = crossSection()/microbarn;
160
161      for (auto& item : _hinc) {
162        if (item.first.substr(0,2) == "04") {
163          double sf = _c["04_sumw"]->val();
164          if (item.first.substr(3,2) == "15")  sf *= 0.8;
165          if (sf)  scale(item.second, xsec/sf);
166          else     normalize(item.second, 0.0);
167        }
168        else {
169          double sf = _c["06_sumw"]->val();
170          if (item.first.substr(3,2) == "15")  sf *= 0.8;
171          if (sf)  scale(item.second, xsec/sf);
172          else     normalize(item.second, 0.0);
173        }
174      }
175
176      for (auto& item : _h) {
177        const string suff = item.first.substr(item.first.size()-14);
178        const double sf = _c["numjets"+suff]->val();
179        if (sf)  scale(item.second, 1.0/sf);
180        else     normalize(item.second, 0.0);
181      }
182
183    }
184
185    /// @}
186
187
188  private:
189
190    /// @name Histograms
191    /// @{
192
193    map<string,CounterPtr> _c;
194    map<string,Histo1DPtr> _h, _hinc;
195
196    /// @}
197
198  };
199
200
201  RIVET_DECLARE_PLUGIN(ATLAS_2011_I919017);
202
203}