1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// Measurement of multi-jet cross sections
 10  class ATLAS_2011_I917599 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I917599);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Projections
 24      const FinalState fs;
 25      FastJets j4(fs, JetAlg::ANTIKT, 0.4);
 26      j4.useInvisibles();
 27      declare(j4, "AntiKtJets04");
 28      FastJets j6(fs, JetAlg::ANTIKT, 0.6);
 29      j6.useInvisibles();
 30      declare(j6, "AntiKtJets06");
 31
 32      // Persistent histograms
 33      book(_h_jet_multi_inclusive, 1, 1, 1);
 34      book(_h_jet_multi_ratio, 2, 1, 1);
 35      _h_jet_pT.resize(4);
 36      book(_h_jet_pT[0] ,3, 1, 1);
 37      book(_h_jet_pT[1] ,4, 1, 1);
 38      book(_h_jet_pT[2] ,5, 1, 1);
 39      book(_h_jet_pT[3] ,6, 1, 1);
 40      book(_h_HT_2 ,7, 1, 1);
 41      book(_h_HT_3 ,8, 1, 1);
 42      book(_h_HT_4 ,9, 1, 1);
 43      //
 44      book(_h_pTlead_R06_60_ratio, 10, 1, 1);
 45      book(_h_pTlead_R06_80_ratio, 11, 1, 1);
 46      book(_h_pTlead_R06_110_ratio, 12, 1, 1);
 47      book(_h_pTlead_R04_60_ratio, 13, 1, 1);
 48      book(_h_pTlead_R04_80_ratio, 14, 1, 1);
 49      book(_h_pTlead_R04_110_ratio, 15, 1, 1);
 50      book(_h_HT2_R06_ratio, 16, 1, 1);
 51      book(_h_HT2_R04_ratio, 17, 1, 1);
 52
 53      // Temporary histograms to be divided for the dsigma3/dsigma2 ratios
 54      book(_h_tmp_pTlead_R06_60_2 , "_pTlead_R06_60_2",  refData(10, 1, 1));
 55      book(_h_tmp_pTlead_R06_80_2 , "_pTlead_R06_80_2",  refData(11, 1, 1));
 56      book(_h_tmp_pTlead_R06_110_2, "_pTlead_R06_110_2", refData(12, 1, 1));
 57      book(_h_tmp_pTlead_R06_60_3 , "_pTlead_R06_60_3",  refData(10, 1, 1));
 58      book(_h_tmp_pTlead_R06_80_3 , "_pTlead_R06_80_3",  refData(11, 1, 1));
 59      book(_h_tmp_pTlead_R06_110_3, "_pTlead_R06_110_3", refData(12, 1, 1));
 60      //
 61      book(_h_tmp_pTlead_R04_60_2 , "_pTlead_R04_60_2",  refData(13, 1, 1));
 62      book(_h_tmp_pTlead_R04_80_2 , "_pTlead_R04_80_2",  refData(14, 1, 1));
 63      book(_h_tmp_pTlead_R04_110_2, "_pTlead_R04_110_2", refData(15, 1, 1));
 64      book(_h_tmp_pTlead_R04_60_3 , "_pTlead_R04_60_3",  refData(13, 1, 1));
 65      book(_h_tmp_pTlead_R04_80_3 , "_pTlead_R04_80_3",  refData(14, 1, 1));
 66      book(_h_tmp_pTlead_R04_110_3, "_pTlead_R04_110_3", refData(15, 1, 1));
 67      //
 68      book(_h_tmp_HT2_R06_2, "_HT2_R06_2", refData(16, 1, 1));
 69      book(_h_tmp_HT2_R06_3, "_HT2_R06_3", refData(16, 1, 1));
 70      book(_h_tmp_HT2_R04_2, "_HT2_R04_2", refData(17, 1, 1));
 71      book(_h_tmp_HT2_R04_3, "_HT2_R04_3", refData(17, 1, 1));
 72    }
 73
 74
 75    /// Perform the per-event analysis
 76    void analyze(const Event& event) {
 77
 78      if (_sedges.empty())  _sedges = _h_jet_multi_inclusive->xEdges();
 79
 80      vector<FourMomentum> jets04;
 81      for (const Jet& jet : apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 60*GeV && Cuts::abseta < 2.8)) {
 82          jets04.push_back(jet.momentum());
 83      }
 84
 85      if (jets04.size() > 1 && jets04[0].pT() > 80.0*GeV) {
 86        for (size_t i = 2; i <= jets04.size(); ++i) {
 87          _h_jet_multi_inclusive->fill(discEdge(i));
 88        }
 89
 90        double HT = 0.0;
 91        for (size_t i = 0; i < jets04.size(); ++i) {
 92          if (i < _h_jet_pT.size()) _h_jet_pT[i]->fill(jets04[i].pT());
 93          HT += jets04[i].pT();
 94        }
 95
 96        if (jets04.size() >= 2) _h_HT_2->fill(HT);
 97        if (jets04.size() >= 3) _h_HT_3->fill(HT);
 98        if (jets04.size() >= 4) _h_HT_4->fill(HT);
 99
100        double pT1(jets04[0].pT()), pT2(jets04[1].pT());
101        double HT2 = pT1 + pT2;
102        if (jets04.size() >= 2) {
103          _h_tmp_HT2_R04_2->fill(HT2);
104          _h_tmp_pTlead_R04_60_2->fill(pT1);
105          if (pT2 > 80.0*GeV) _h_tmp_pTlead_R04_80_2->fill(pT1);
106          if (pT2 > 110.0*GeV) _h_tmp_pTlead_R04_110_2->fill(pT1);
107        }
108        if (jets04.size() >= 3) {
109          double pT3(jets04[2].pT());
110          _h_tmp_HT2_R04_3->fill(HT2);
111          _h_tmp_pTlead_R04_60_3->fill(pT1);
112          if (pT3 > 80.0*GeV) _h_tmp_pTlead_R04_80_3->fill(pT1);
113          if (pT3 > 110.0*GeV) _h_tmp_pTlead_R04_110_3->fill(pT1);
114        }
115      }
116
117      /// @todo It'd be better to avoid duplicating 95% of the code!
118      vector<FourMomentum> jets06;
119      for (const Jet& jet : apply<FastJets>(event, "AntiKtJets06").jetsByPt(Cuts::pT > 60*GeV && Cuts::abseta < 2.8)) {
120          jets06.push_back(jet.momentum());
121      }
122      if (jets06.size() > 1 && jets06[0].pT() > 80.0*GeV) {
123        double pT1(jets06[0].pT()), pT2(jets06[1].pT());
124        double HT2 = pT1 + pT2;
125        if (jets06.size() >= 2) {
126          _h_tmp_HT2_R06_2->fill(HT2);
127          _h_tmp_pTlead_R06_60_2->fill(pT1);
128          if (pT2 > 80.0*GeV) _h_tmp_pTlead_R06_80_2->fill(pT1);
129          if (pT2 > 110.0*GeV) _h_tmp_pTlead_R06_110_2->fill(pT1);
130        }
131        if (jets06.size() >= 3) {
132          double pT3(jets06[2].pT());
133          _h_tmp_HT2_R06_3->fill(HT2);
134          _h_tmp_pTlead_R06_60_3->fill(pT1);
135          if (pT3 > 80.0*GeV) _h_tmp_pTlead_R06_80_3->fill(pT1);
136          if (pT3 > 110.0*GeV) _h_tmp_pTlead_R06_110_3->fill(pT1);
137        }
138      }
139
140    }
141
142
143    /// Normalise histograms etc., after the run
144    void finalize() {
145
146      // Normalize std histos
147      scale(_h_jet_multi_inclusive, crossSectionPerEvent());
148      scale(_h_jet_pT[0], crossSectionPerEvent());
149      scale(_h_jet_pT[1], crossSectionPerEvent());
150      scale(_h_jet_pT[2], crossSectionPerEvent());
151      scale(_h_jet_pT[3], crossSectionPerEvent());
152      scale(_h_HT_2, crossSectionPerEvent());
153      scale(_h_HT_3, crossSectionPerEvent());
154      scale(_h_HT_4, crossSectionPerEvent());
155
156      scale(_h_tmp_pTlead_R06_60_2 , crossSectionPerEvent());
157      scale(_h_tmp_pTlead_R06_80_2 , crossSectionPerEvent());
158      scale(_h_tmp_pTlead_R06_110_2, crossSectionPerEvent());
159      scale(_h_tmp_pTlead_R06_60_3 , crossSectionPerEvent());
160      scale(_h_tmp_pTlead_R06_80_3 , crossSectionPerEvent());
161      scale(_h_tmp_pTlead_R06_110_3, crossSectionPerEvent());
162
163      scale(_h_tmp_pTlead_R04_60_2 , crossSectionPerEvent());
164      scale(_h_tmp_pTlead_R04_80_2 , crossSectionPerEvent());
165      scale(_h_tmp_pTlead_R04_110_2, crossSectionPerEvent());
166      scale(_h_tmp_pTlead_R04_60_3 , crossSectionPerEvent());
167      scale(_h_tmp_pTlead_R04_80_3 , crossSectionPerEvent());
168      scale(_h_tmp_pTlead_R04_110_3, crossSectionPerEvent());
169
170      scale(_h_tmp_HT2_R06_2, crossSectionPerEvent());
171      scale(_h_tmp_HT2_R06_3, crossSectionPerEvent());
172      scale(_h_tmp_HT2_R04_2, crossSectionPerEvent());
173      scale(_h_tmp_HT2_R04_3, crossSectionPerEvent());
174
175      // Fill inclusive jet multiplicity ratio
176      for (auto& b : _h_jet_multi_ratio->bins()) {
177        size_t idx = b.index();
178        if (_h_jet_multi_inclusive->bin(idx).sumW()) {
179          const double val = _h_jet_multi_inclusive->bin(idx+1).sumW() / _h_jet_multi_inclusive->bin(idx).sumW();
180          // @todo Shouldn't these be added in quadrature??
181          const double err = ( _h_jet_multi_inclusive->bin(idx+1).relErrW() + _h_jet_multi_inclusive->bin(idx).relErrW() ) * val;
182          b.set(val, err);
183        }
184      }
185
186      /// Create ratio histograms
187      divide(_h_tmp_pTlead_R06_60_3,_h_tmp_pTlead_R06_60_2, _h_pTlead_R06_60_ratio);
188      divide(_h_tmp_pTlead_R06_80_3,_h_tmp_pTlead_R06_80_2, _h_pTlead_R06_80_ratio);
189      divide(_h_tmp_pTlead_R06_110_3,_h_tmp_pTlead_R06_110_2, _h_pTlead_R06_110_ratio);
190      divide(_h_tmp_pTlead_R04_60_3,_h_tmp_pTlead_R04_60_2, _h_pTlead_R04_60_ratio);
191      divide(_h_tmp_pTlead_R04_80_3,_h_tmp_pTlead_R04_80_2, _h_pTlead_R04_80_ratio);
192      divide(_h_tmp_pTlead_R04_110_3,_h_tmp_pTlead_R04_110_2, _h_pTlead_R04_110_ratio);
193      divide(_h_tmp_HT2_R06_3,_h_tmp_HT2_R06_2, _h_HT2_R06_ratio);
194      divide(_h_tmp_HT2_R04_3,_h_tmp_HT2_R04_2, _h_HT2_R04_ratio);
195    }
196
197    string discEdge(size_t m) const {
198      size_t idx = m - 2;
199      if (idx < _sedges.size())  return _sedges[idx];
200      return "OTHER";
201    }
202
203    /// @}
204
205
206  private:
207
208    /// @name Histograms
209    /// @{
210    BinnedHistoPtr<string> _h_jet_multi_inclusive;
211    BinnedEstimatePtr<string> _h_jet_multi_ratio;
212    vector<string> _sedges;
213
214    vector<Histo1DPtr> _h_jet_pT;
215    Histo1DPtr _h_HT_2;
216    Histo1DPtr _h_HT_3;
217    Histo1DPtr _h_HT_4;
218    /// @}
219
220    /// @name Ratio histograms
221    /// @{
222    Estimate1DPtr _h_pTlead_R06_60_ratio, _h_pTlead_R06_80_ratio, _h_pTlead_R06_110_ratio;
223    Estimate1DPtr _h_pTlead_R04_60_ratio, _h_pTlead_R04_80_ratio, _h_pTlead_R04_110_ratio;
224    Estimate1DPtr _h_HT2_R06_ratio, _h_HT2_R04_ratio;
225    /// @}
226
227    /// @name Temporary histograms to be divided for the dsigma3/dsigma2 ratios
228    /// @{
229    Histo1DPtr _h_tmp_pTlead_R06_60_2, _h_tmp_pTlead_R06_80_2, _h_tmp_pTlead_R06_110_2;
230    Histo1DPtr _h_tmp_pTlead_R06_60_3, _h_tmp_pTlead_R06_80_3, _h_tmp_pTlead_R06_110_3;
231    Histo1DPtr _h_tmp_pTlead_R04_60_2, _h_tmp_pTlead_R04_80_2, _h_tmp_pTlead_R04_110_2;
232    Histo1DPtr _h_tmp_pTlead_R04_60_3, _h_tmp_pTlead_R04_80_3, _h_tmp_pTlead_R04_110_3;
233    Histo1DPtr _h_tmp_HT2_R06_2, _h_tmp_HT2_R06_3, _h_tmp_HT2_R04_2, _h_tmp_HT2_R04_3;
234    /// @}
235
236  };
237
238
239  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I917599, ATLAS_2011_S9128077);
240
241}