1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Tools/Logging.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Tools/BinnedHistogram.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief ATLAS azimuthal decorrelation with jet veto analysis
 11  /// @author James Robinson <james.robinson@cern.ch>
 12  class ATLAS_2014_I1307243 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2014_I1307243);
 17
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      _dy_max = 8; _years = { 2010, 2011 };
 23      _Qnoughts = { 20., 30., 40., 50., 60., 70., 80., 90., 100. };
 24      /// Initialise and register projections here
 25      FastJets fastJets(FinalState(), FastJets::ANTIKT, 0.6, JetAlg::Muons::ALL, JetAlg::Invisibles::ALL);
 26      declare(fastJets, "AntiKt6JetsWithInvisibles");
 27
 28
 29      /// Book histograms
 30      for (const string cat : { "inclusive", "gap" }) {
 31        const size_t offset = (cat == "gap") ? 1 : 0;
 32
 33        // Temporary inclusive and gap histograms
 34        book(_aux_dy[cat], "_" + cat + "_dy", refData(1, 1, 1));
 35        book(_aux_pTbar[cat], "_" + cat + "_pTbar", refData(2, 1, 1));
 36
 37        book(_h_C2C1_dy[cat],    7 + 4 * offset, 1, 1, false);
 38        book(_h_C2C1_pTbar[cat], 8 + 4 * offset, 1, 1, false);
 39
 40        // Azimuthal moment histograms
 41        book(_p_cosDeltaPhi_dy[cat],        5 + 4 * offset, 1, 1);
 42        book(_p_cosDeltaPhi_pTbar[cat],     6 + 4 * offset, 1, 1);
 43        book(_p_cosTwoDeltaPhi_dy[cat],    37 + 2 * offset, 1, 1);
 44        book(_p_cosTwoDeltaPhi_pTbar[cat], 38 + 2 * offset, 1, 1);
 45
 46        // Gap fraction vs. Q0 and cross-section in dy slices
 47        _s_gapFrac_Q0.resize(_dy_max);
 48        for (size_t dyLow = 0; dyLow < _dy_max; ++dyLow ) {
 49          const string hname("_" + cat + "_dySlice_" + toString(dyLow) + "_" + toString(dyLow+1) + "_Q0");
 50          { Histo1DPtr tmp; _aux_Q0_dySlices[cat].add(dyLow, dyLow+1, book(tmp, hname, refData(29+dyLow, 1, 1))); }
 51          { Histo1DPtr tmp; _h_dphi_dySlices[cat].add(dyLow, dyLow+1, book(tmp, 13+(_dy_max*offset)+dyLow, 1, 1)); }
 52          if (!offset)  book(_s_gapFrac_Q0[dyLow], 29 + dyLow, 1, 1); // only book once
 53        }
 54
 55      }
 56
 57      // Number of jets in rapidity interval
 58      book(_s_gapFrac_dy,     1, 1, 1);
 59      book(_s_gapFrac_pTbar,  2, 1, 1);
 60      book(_p_nGapJets_dy,    3, 1, 1);
 61      book(_p_nGapJets_pTbar, 4, 1, 1);
 62    }
 63
 64
 65    /// Perform the per-event analysis
 66    void analyze(const Event& event) {
 67
 68      for (size_t reg = 0; reg < 2; ++reg) {
 69
 70        // Retrieve all anti-kt R=0.6 jets 
 71        const double maxRap = reg? 2.4 : 4.4;
 72        const Jets& akt6Jets = apply<JetAlg>(event, "AntiKt6JetsWithInvisibles").jetsByPt(Cuts::absrap < maxRap);
 73        // If there are fewer than 2 jets then bail
 74        if ( akt6Jets.size() < 2 ) { vetoEvent; }
 75
 76        // Require jets to be above {60, 50} GeV
 77        if ( akt6Jets[0].pT() < 60*GeV || akt6Jets[1].pT() < 50*GeV ) { vetoEvent; }
 78
 79        // Identify gap boundaries
 80        const double yMin = std::min(akt6Jets[0].rap(), akt6Jets[1].rap());
 81        const double yMax = std::max(akt6Jets[0].rap(), akt6Jets[1].rap());
 82
 83        // Determine azimuthal decorrelation quantities
 84        const double dy = yMax - yMin;
 85        const double dphi = mapAngle0ToPi(deltaPhi(akt6Jets[0], akt6Jets[1]));
 86        const double pTbar = 0.5*(akt6Jets[0].pT() + akt6Jets[1].pT())/GeV;
 87
 88        // Impose minimum dy for the 2011 phase space
 89        if ( _years[reg] == 2011 && dy < 1.0 ) { vetoEvent; }
 90
 91        // Determine gap quantities
 92        size_t nGapJets = 0;
 93        double maxGapQ0 = 0.;
 94        const double vetoScale = reg? 30*GeV : 20*GeV;
 95        for (const Jet& jet : akt6Jets) {
 96          if (!inRange(jet.rap(), yMin, yMax, OPEN, OPEN))  continue;
 97          const double pT = jet.pT()/GeV;
 98          if (pT > vetoScale) { ++nGapJets; }
 99          if (pT > maxGapQ0) { maxGapQ0 = pT; }
100        }
101
102        // Fill histograms
103        fillHists(_years[reg], nGapJets, { dy, pTbar, dphi, maxGapQ0 });
104      }
105      return;
106    }
107
108    void fillHists(const size_t region, const size_t nGapJets, const vector<double>& vars) {
109      assert(vars.size() == 4);
110      const double dy = vars[0];
111      const double pTbar = vars[1];
112      const double dphi = vars[2];
113      const double maxGapQ0 = vars[3];
114      // Determine gap category
115      vector<string> categories = {"inclusive"};
116      if (!nGapJets) { categories += string("gap"); }
117
118      // Fill histograms relevant for comparison with 2010 data
119      if (region == _years[0]) {
120        // Fill inclusive and gap histograms
121        for (const string& cat : categories) {
122          _aux_dy[cat]->fill(dy);
123          _h_dphi_dySlices[cat].fill(dy, dphi/M_PI);
124          _p_cosDeltaPhi_dy[cat]->fill(dy, cos(M_PI - dphi));
125          _p_cosTwoDeltaPhi_dy[cat]->fill(dy, cos(2*dphi));
126        }
127        // Fill profiled nGapJets
128        _p_nGapJets_dy->fill(dy, nGapJets);
129        // Fill Q0 histograms - can fill multiple points per event
130        for (const double Q0 :  _Qnoughts) {
131          _aux_Q0_dySlices["inclusive"].fill(dy, Q0);
132          if (maxGapQ0 <= Q0) { _aux_Q0_dySlices["gap"].fill(dy, Q0); }
133        }
134
135      // Fill histograms relevant for comparison with 2011 data
136      } 
137      else if (region == _years[1]) {
138        // Fill inclusive and gap histograms
139        for (const string& cat : categories) {
140          _aux_pTbar[cat]->fill(pTbar);
141          _p_cosDeltaPhi_pTbar[cat]->fill(pTbar, cos(M_PI - dphi));
142          _p_cosTwoDeltaPhi_pTbar[cat]->fill(pTbar, cos(2*dphi));
143        }
144        // Fill profiled nGapJets
145        _p_nGapJets_pTbar->fill(pTbar, nGapJets);
146      }
147      return;
148    }
149
150    /// Normalise histograms etc., after the run
151    void finalize() {
152      // Normalise cross-section plots to correct cross-section
153      const double ySpan = 1.0; // all dy spans are 1
154      const double sf = crossSection() / picobarn / sumOfWeights();
155      for (const string cat : { "inclusive", "gap" }) {
156        _h_dphi_dySlices[cat].scale(sf/ySpan/M_PI, this);
157        // Create C2/C1 scatter from profiles
158        divide(_p_cosTwoDeltaPhi_dy[cat],    _p_cosDeltaPhi_dy[cat],    _h_C2C1_dy[cat]);
159        divide(_p_cosTwoDeltaPhi_pTbar[cat], _p_cosDeltaPhi_pTbar[cat], _h_C2C1_pTbar[cat]);
160      }
161
162      // Fill simple gap fractions
163      efficiency(_aux_dy["gap"], _aux_dy["inclusive"], _s_gapFrac_dy);
164      efficiency(_aux_pTbar["gap"], _aux_pTbar["inclusive"], _s_gapFrac_pTbar);
165
166      // Register and fill Q0 gap fractions
167      for (size_t dyLow = 0; dyLow < _dy_max; ++dyLow) {
168        efficiency(_aux_Q0_dySlices["gap"].histos()[dyLow], _aux_Q0_dySlices["inclusive"].histos()[dyLow], _s_gapFrac_Q0[dyLow]);
169      }
170    }
171
172
173  private:
174
175    /// Member variables
176    vector<size_t> _years;
177    vector<double> _Qnoughts;
178
179    size_t _dy_max;
180
181    /// auxiliary histograms for gap fractions
182    map<string, Histo1DPtr> _aux_dy;
183    map<string, Histo1DPtr> _aux_pTbar;
184    map<string, BinnedHistogram> _aux_Q0_dySlices;
185
186    // Gap fractions
187    Scatter2DPtr _s_gapFrac_dy, _s_gapFrac_pTbar;
188    vector<Scatter2DPtr> _s_gapFrac_Q0;
189
190    // Number of jets in rapidity interval
191    Profile1DPtr _p_nGapJets_dy;
192    Profile1DPtr _p_nGapJets_pTbar;
193
194    // Azimuthal moment histograms
195    map<string, Profile1DPtr> _p_cosDeltaPhi_dy;
196    map<string, Profile1DPtr> _p_cosDeltaPhi_pTbar;
197    map<string, Profile1DPtr> _p_cosTwoDeltaPhi_dy;
198    map<string, Profile1DPtr> _p_cosTwoDeltaPhi_pTbar;
199    map<string, Scatter2DPtr> _h_C2C1_dy;
200    map<string, Scatter2DPtr> _h_C2C1_pTbar;
201
202    // Cross-section vs. deltaPhi in deltaY slices
203    map<string, BinnedHistogram> _h_dphi_dySlices;
204  };
205
206  // The hook for the plugin system
207  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1307243);
208
209}