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