rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2011_S9126244

Measurement of dijet production with a veto on additional central jet activity
Experiment: ATLAS (LHC)
Inspire ID: 917526
Status: VALIDATED
Authors:
  • Graham Jones
References: Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • Require QCD interactions at 7TeV. A substantial number of events are required to populate the large rapidity seperation region.

A measurement of the jet activity in rapidity intervals bounded by a dijet system. The fraction of events passing a veto requirement are shown as a function of both the rapidity interval size and the average transverse momentum of the dijet system. The average number of jets above the veto threshold are also shown as a function of the same variables. There are two possible selection criteria applied to data. Either the two highest transverse momentum jets or the jets most forward and backward in rapidity are taken to define the dijet system, where the veto threhsold is 20 GeV. Additionally for the latter selection an alternative veto transverse momentum threshold which is equal to the average transverse momentum is applied. Jet selections are based on the anti-$k_t$ algorithm with $R=0.6$, $p_\perp > 20$ GeV and $|y_\text{jet}| < 4.4$.

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