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  struct ATLAS_2011_I917526_Plots {
 10
 11    int selectionType; ///< The HepData y-axis code
 12    string intermediateHistName;
 13
 14    // Gap fraction vs DeltaY plot setup
 15    int _gapFractionDeltaYHistIndex;
 16    vector<double> _gapFractionDeltaYSlices;
 17    Histo1DGroupPtr _h_gapVsDeltaYVeto;
 18    Histo1DGroupPtr _h_gapVsDeltaYInc;
 19    vector<Estimate1DPtr> _ratio_DeltaY;
 20
 21    // Gap fraction vs ptBar plot setup
 22    int _gapFractionPtBarHistIndex;
 23    vector<double> _gapFractionPtBarSlices;
 24    Histo1DGroupPtr _h_gapVsPtBarVeto;
 25    Histo1DGroupPtr _h_gapVsPtBarInc;
 26    vector<Estimate1DPtr> _ratio_PtBar;
 27
 28    // Gap fraction vs Q0 plot setup
 29    int _gapFractionQ0HistIndex;
 30    vector<double> _gapFractionQ0SlicesPtBar;
 31    vector<double> _gapFractionQ0SlicesDeltaY;
 32    vector<Histo1DPtr> _h_vetoPt;
 33    vector<Estimate1DPtr> _d_vetoPtGapFraction;
 34    vector<double> _vetoPtTotalSum; ///< @todo Can this just be replaced with _h_vetoPt.integral()?
 35
 36    // Average njet vs DeltaY setup
 37    int _avgNJetDeltaYHistIndex;
 38    vector<double> _avgNJetDeltaYSlices;
 39    vector<Profile1DPtr> _p_avgJetVsDeltaY;
 40
 41    // Average njet vs PptBar setup
 42    int _avgNJetPtBarHistIndex;
 43    vector<double> _avgNJetPtBarSlices;
 44    vector<Profile1DPtr> _p_avgJetVsPtBar;
 45  };
 46
 47
 48
 49  /// ATLAS dijet production with central jet veto
 50  ///
 51  /// @todo Make sure that temp histos are removed
 52  class ATLAS_2011_I917526 : public Analysis {
 53  public:
 54
 55    /// Constructor
 56    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I917526);
 57
 58
 59    /// Book histograms and initialise projections before the run
 60    void init() {
 61
 62      // Initialize the lone projection required
 63      declare(FastJets(FinalState(), JetAlg::ANTIKT, 0.6), "AntiKtJets06");
 64
 65      // Initialize plots for each selection type
 66      _selectionPlots[0].intermediateHistName = "highestPt";
 67      _selectionPlots[0].selectionType = 1;
 68      _selectionPlots[0]._gapFractionDeltaYHistIndex = 6;
 69      _selectionPlots[0]._gapFractionPtBarHistIndex = 1;
 70      _selectionPlots[0]._gapFractionQ0HistIndex = 13;
 71      _selectionPlots[0]._avgNJetDeltaYHistIndex = 37;
 72      _selectionPlots[0]._avgNJetPtBarHistIndex = 26;
 73      _selectionPlots[0]._gapFractionDeltaYSlices = {{ 70.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0 }};
 74      _selectionPlots[0]._gapFractionPtBarSlices = {{ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }};
 75      _selectionPlots[0]._gapFractionQ0SlicesPtBar = {{ 70.0, 90.0, 120.0, 150.0, 210.0, 240.0 }};
 76      _selectionPlots[0]._gapFractionQ0SlicesDeltaY = {{ 2.0, 3.0, 4.0, 5.0 }};
 77      _selectionPlots[0]._avgNJetPtBarSlices = {{ 1.0, 2.0, 3.0, 4.0, 5.0 }};
 78      _selectionPlots[0]._avgNJetDeltaYSlices = {{ 70.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0 }};
 79      initializePlots(_selectionPlots[0]);
 80
 81      _selectionPlots[1].intermediateHistName = "forwardBackward";
 82      _selectionPlots[1].selectionType = 2;
 83      _selectionPlots[1]._gapFractionDeltaYHistIndex = 6;
 84      _selectionPlots[1]._gapFractionPtBarHistIndex = 1;
 85      _selectionPlots[1]._gapFractionQ0HistIndex = 13;
 86      _selectionPlots[1]._avgNJetDeltaYHistIndex = 37;
 87      _selectionPlots[1]._avgNJetPtBarHistIndex = 26;
 88      _selectionPlots[1]._gapFractionDeltaYSlices = {{ 70.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0 }};
 89      _selectionPlots[1]._gapFractionPtBarSlices = {{ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }};
 90      _selectionPlots[1]._gapFractionQ0SlicesPtBar = {{ 70.0, 90.0, 120.0, 150.0, 210.0, 240.0 }};
 91      _selectionPlots[1]._gapFractionQ0SlicesDeltaY = {{ 2.0, 3.0, 4.0, 5.0 }};
 92      _selectionPlots[1]._avgNJetPtBarSlices = {{ 1.0, 2.0, 3.0, 4.0, 5.0 }};
 93      _selectionPlots[1]._avgNJetDeltaYSlices = {{ 70.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0 }};
 94      initializePlots(_selectionPlots[1]);
 95
 96      _selectionPlots[2].intermediateHistName = "forwardBackward_PtBarVeto";
 97      _selectionPlots[2].selectionType = 1;
 98      _selectionPlots[2]._gapFractionDeltaYHistIndex = 19;
 99      _selectionPlots[2]._avgNJetDeltaYHistIndex = 30;
100      _selectionPlots[2]._gapFractionDeltaYSlices = {{ 70.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0 }};
101      _selectionPlots[2]._avgNJetDeltaYSlices = {{ 70.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0 }};
102      initializePlots(_selectionPlots[2]);
103    }
104
105
106    void initializePlots(ATLAS_2011_I917526_Plots& plots) {
107
108      // Gap fraction vs DeltaY
109      if (!plots._gapFractionDeltaYSlices.empty()) {
110        book(plots._h_gapVsDeltaYVeto, plots._gapFractionDeltaYSlices);
111        book(plots._h_gapVsDeltaYInc,  plots._gapFractionDeltaYSlices);
112        for (size_t x = 0; x < plots._h_gapVsDeltaYInc->numBins(); ++x) {
113          const string vetoHistName = "TMP/gapDeltaYVeto_" + plots.intermediateHistName + "_" + to_str(x);
114          const string inclusiveHistName = "TMP/gapDeltaYInclusive_" + plots.intermediateHistName + "_" + to_str(x);
115          const auto& ref = refData(plots._gapFractionDeltaYHistIndex+x, 1,plots.selectionType);
116          book(plots._h_gapVsDeltaYVeto->bin(x+1), vetoHistName, ref);
117          book(plots._h_gapVsDeltaYInc->bin(x+1), inclusiveHistName, ref);
118          plots._ratio_DeltaY.push_back(Estimate1DPtr());
119          book(plots._ratio_DeltaY[x], plots._gapFractionDeltaYHistIndex+x, 1, plots.selectionType);
120        }
121      }
122
123      // Average njet vs DeltaY
124      if (!plots._avgNJetDeltaYSlices.empty()) {
125        for (size_t x = 0; x < plots._avgNJetDeltaYSlices.size()-1; x++) {
126          Profile1DPtr tmp;
127          plots._p_avgJetVsDeltaY += book(tmp, plots._avgNJetDeltaYHistIndex+x, 1, plots.selectionType);
128        }
129      }
130
131      // Gap fraction vs PtBar
132      if (!plots._gapFractionPtBarSlices.empty()) {
133
134        book(plots._h_gapVsPtBarVeto, plots._gapFractionPtBarSlices);
135        book(plots._h_gapVsPtBarInc,  plots._gapFractionPtBarSlices);
136        for (size_t x = 0; x < plots._h_gapVsPtBarInc->numBins(); ++x) {
137          const string vetoHistName = "TMP/gapPtBarVeto_" + plots.intermediateHistName + "_" + to_str(x);
138          const string inclusiveHistName = "TMP/gapPtBarInclusive_" + plots.intermediateHistName + "_" + to_str(x);
139          const auto& ref = refData(plots._gapFractionPtBarHistIndex+x, 1, plots.selectionType);
140          book(plots._h_gapVsPtBarVeto->bin(x+1), vetoHistName, ref);
141          book(plots._h_gapVsPtBarInc->bin(x+1), inclusiveHistName, ref);
142          plots._ratio_PtBar.push_back(Estimate1DPtr());
143          book(plots._ratio_PtBar[x], plots._gapFractionPtBarHistIndex+x, 1, plots.selectionType);
144        }
145      }
146
147      // Average njet vs PtBar
148      if (!plots._avgNJetPtBarSlices.empty()) {
149        for (size_t x=0; x<plots._avgNJetPtBarSlices.size()-1; x++) {
150          Profile1DPtr tmp;
151          plots._p_avgJetVsPtBar += book(tmp, plots._avgNJetPtBarHistIndex+x, 1, plots.selectionType);
152        }
153      }
154
155      // Gap fraction vs Q0
156      int q0PlotCount = 0;
157      for (size_t x = 0; x < plots._gapFractionQ0SlicesPtBar.size()/2; x++) {
158        for (size_t y = 0; y < plots._gapFractionQ0SlicesDeltaY.size()/2; y++) {
159          const string vetoPtHistName = "TMP/vetoPt_" + plots.intermediateHistName + "_" + to_str(q0PlotCount);
160          Histo1DPtr tmp1;
161          plots._h_vetoPt += book(tmp1,vetoPtHistName, refData(plots._gapFractionQ0HistIndex + q0PlotCount, 1, plots.selectionType));
162          Estimate1DPtr tmp2;
163          plots._d_vetoPtGapFraction += book(tmp2,plots._gapFractionQ0HistIndex + q0PlotCount, 1, plots.selectionType);
164          plots._vetoPtTotalSum += 0.0; ///< @todo Can this just be replaced with _h_vetoPt.integral()?
165          ++q0PlotCount;
166        }
167      }
168    }
169
170
171    /// Perform the per-event analysis
172    void analyze(const Event& event) {
173
174      // Get minimal list of jets needed to be considered
175      double minimumJetPtBar = 50.0*GeV; // of interval defining jets
176
177      vector<FourMomentum> acceptJets;
178      for (const Jet& jet : apply<FastJets>(event, "AntiKtJets06").jetsByPt(Cuts::pT > 20*GeV && Cuts::absrap < 4.4)) {
179        acceptJets.push_back(jet.momentum());
180      }
181
182      // If we can't form an interval, drop out of the analysis early
183      if (acceptJets.size() < 2) vetoEvent;
184
185      // Analyze leading jet case
186      if (acceptJets[0].pT() + acceptJets[1].pT() > 2*minimumJetPtBar) {
187        analyzeJets(acceptJets, _selectionPlots[0], 20.0*GeV);
188      }
189
190      // Find the most forward-backward jets
191      size_t minRapidityJet = 0, maxRapidityJet = 0;
192      for (size_t j = 1; j < acceptJets.size(); j++) {
193        if (acceptJets[j].rapidity() > acceptJets[maxRapidityJet].rapidity()) maxRapidityJet = j;
194        if (acceptJets[j].rapidity() < acceptJets[minRapidityJet].rapidity()) minRapidityJet = j;
195      }
196
197      // Make a container of jet momenta with the extreme f/b jets at the front
198      vector<FourMomentum> fwdBkwdJets;
199      fwdBkwdJets.push_back(acceptJets[maxRapidityJet]);
200      fwdBkwdJets.push_back(acceptJets[minRapidityJet]);
201      for (size_t j = 0; j < acceptJets.size(); j++) {
202        if (j == minRapidityJet || j == maxRapidityJet) continue;
203        fwdBkwdJets.push_back(acceptJets[j]);
204      }
205
206      if (fwdBkwdJets[0].pT() + fwdBkwdJets[1].pT() > 2*minimumJetPtBar) {
207        // Use most forward/backward jets in rapidity to define the interval
208        analyzeJets(fwdBkwdJets, _selectionPlots[1], 20.0*GeV);
209        // As before but now using PtBar of interval to define veto threshold
210        analyzeJets(fwdBkwdJets, _selectionPlots[2], (fwdBkwdJets[0].pT()+fwdBkwdJets[1].pT())/2.0);
211      }
212    }
213
214
215    /// Fill plots!
216    void analyzeJets(vector<FourMomentum>& jets, ATLAS_2011_I917526_Plots& plots,
217                     double vetoPtThreshold) {
218
219      // Calculate the interval size, ptBar and veto Pt (if any)
220      const double intervalSize = fabs(jets[0].rapidity()-jets[1].rapidity());
221      const double ptBar = (jets[0].pT()+jets[1].pT())/2.0;
222
223      const double minY = min(jets[0].rapidity(), jets[1].rapidity());
224      const double maxY = max(jets[0].rapidity(), jets[1].rapidity());
225
226      double vetoPt = 0.0*GeV;
227      for (size_t j = 2; j < jets.size(); j++) {
228        if (inRange(jets[j].rapidity(), minY, maxY)) vetoPt = max(jets[j].pT(), vetoPt);
229      }
230
231      // Fill the gap fraction vs delta Y histograms
232      if (plots._gapFractionDeltaYSlices.size()) {
233        plots._h_gapVsDeltaYInc->fill(ptBar/GeV, intervalSize);
234        if (vetoPt < vetoPtThreshold) {
235          plots._h_gapVsDeltaYVeto->fill(ptBar/GeV, intervalSize);
236        }
237      }
238
239      // Fill the gap fraction vs ptBar histograms
240      if (plots._gapFractionPtBarSlices.size()) {
241         plots._h_gapVsPtBarInc->fill(intervalSize, ptBar/GeV);
242        if (vetoPt < vetoPtThreshold) {
243          plots._h_gapVsPtBarVeto->fill(intervalSize, ptBar/GeV);
244        }
245      }
246
247      // Count the number of veto jets present
248      int vetoJetsCount = 0;
249      for (size_t j = 2; j < jets.size(); ++j) {
250        if (inRange(jets[j].rapidity(), minY, maxY) && jets[j].pT() > vetoPtThreshold) {
251          vetoJetsCount += 1;
252        }
253      }
254
255      // Fill the avg NJet, deltaY slices
256      if (!plots._avgNJetPtBarSlices.empty()) {
257        for (size_t i = 0; i < plots._avgNJetPtBarSlices.size()-1; ++i) {
258          if (inRange(intervalSize, plots._avgNJetPtBarSlices[i], plots._avgNJetPtBarSlices[i+1])) {
259            plots._p_avgJetVsPtBar[i]->fill(ptBar/GeV, vetoJetsCount);
260          }
261        }
262      }
263
264      // Fill the avg NJet, ptBar slices
265      if (!plots._avgNJetDeltaYSlices.empty()) {
266        for (size_t i = 0; i < plots._avgNJetDeltaYSlices.size()-1; ++i) {
267          if (inRange(ptBar/GeV, plots._avgNJetDeltaYSlices[i], plots._avgNJetDeltaYSlices[i+1])) {
268            plots._p_avgJetVsDeltaY[i]->fill(intervalSize, vetoJetsCount);
269          }
270        }
271      }
272
273      // Fill the veto pt plots
274      int q0PlotCount = 0;
275      for (size_t x = 0; x < plots._gapFractionQ0SlicesPtBar.size()/2; ++x) {
276        for (size_t y = 0; y < plots._gapFractionQ0SlicesDeltaY.size()/2; ++y) {
277          // Check if it should be filled
278          if ( ptBar/GeV < plots._gapFractionQ0SlicesPtBar[x*2] ||
279               ptBar/GeV >= plots._gapFractionQ0SlicesPtBar[x*2+1] ) {
280            q0PlotCount++;
281            continue;
282          }
283
284          if ( intervalSize < plots._gapFractionQ0SlicesDeltaY[y*2] ||
285               intervalSize >= plots._gapFractionQ0SlicesDeltaY[y*2+1] ) {
286            q0PlotCount++;
287            continue;
288          }
289
290          plots._h_vetoPt[q0PlotCount]->fill(vetoPt);
291          plots._vetoPtTotalSum[q0PlotCount] += 1.0;
292
293          q0PlotCount++;
294        }
295      }
296    }
297
298
299    /// Derive final distributions for each selection
300    void finalize() {
301
302      for (const ATLAS_2011_I917526_Plots& plots : _selectionPlots) {
303
304        /// @todo Clean up temp histos -- requires restructuring the temp histo struct
305
306        if (plots._gapFractionDeltaYSlices.size()) {
307          for (size_t x = 0; x < plots._h_gapVsDeltaYVeto->numBins(); ++x) {
308            divide(plots._h_gapVsDeltaYVeto->bin(x+1), plots._h_gapVsDeltaYInc->bin(x+1), plots._ratio_DeltaY[x]);
309          }
310        }
311
312        if (plots._gapFractionPtBarSlices.size()) {
313          for (size_t x = 0; x < plots._h_gapVsPtBarVeto->numBins(); ++x) {
314            divide(plots._h_gapVsPtBarVeto->bin(x+1), plots._h_gapVsPtBarInc->bin(x+1), plots._ratio_PtBar[x]);
315          }
316        }
317
318        for (size_t h = 0; h < plots._d_vetoPtGapFraction.size(); ++h) {
319          finalizeQ0GapFraction(plots._vetoPtTotalSum[h], plots._d_vetoPtGapFraction[h], plots._h_vetoPt[h]);
320        }
321      }
322    }
323
324
325    /// Convert the differential histograms to an integral histo and assign binomial errors as a efficiency
326    /// @todo Should be convertible to a YODA ~one-liner using toIntegralEfficiencyHisto
327    void finalizeQ0GapFraction(double totalWeightSum, Estimate1DPtr gapFractionDP, const Histo1DPtr& vetoPtHist) {
328      for (const auto& b : vetoPtHist->bins()) {
329        const double vetoPtWeightSum = vetoPtHist->integralTo(b.index()-1); ///< Integral (with underflow) up to but not including bin i
330        // Calculate the efficiency & binomial uncertainty
331        const double eff = (totalWeightSum != 0) ? vetoPtWeightSum/totalWeightSum : 0;
332        const double effErr = (totalWeightSum != 0) ? sqrt( eff*(1.0-eff)/totalWeightSum ) : 0;
333        gapFractionDP->bin(b.index()).set(eff, effErr);
334      }
335    }
336
337
338  private:
339
340    // Struct containing the complete set of plots, times 3 for the different selections
341    ATLAS_2011_I917526_Plots _selectionPlots[3];
342
343  };
344
345
346  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I917526, ATLAS_2011_S9126244);
347
348}