Rivet analyses referenceATLAS_2011_S9126244Measurement of dijet production with a veto on additional central jet activityExperiment: ATLAS (LHC) Inspire ID: 917526 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
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}
|