Rivet analyses referenceATLAS_2011_I945498$Z$+jets in $pp$ at 7 TeVExperiment: ATLAS (LHC) Inspire ID: 945498 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
Production of jets in association with a $Z/\gamma^*$ boson in proton--proton collisions at $\sqrt{s} = 7$ TeV with the ATLAS detector. The analysis includes the full 2010 data set, collected with a low rate of multiple proton--proton collisions in the accelerator, corresponding to an integrated luminosity of 36 pb$^{-1}$. Inclusive jet cross sections in $Z/\gamma^*$ events, with $Z/\gamma^*$ decaying into electron or muon pairs, are measured for jets with transverse momentum $p_T > 30$ GeV and jet rapidity $|y| < 4.4$. Source code: ATLAS_2011_I945498.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/DileptonFinder.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/FinalState.hh"
6#include "Rivet/Projections/VetoedFinalState.hh"
7#include "Rivet/Projections/IdentifiedFinalState.hh"
8#include "Rivet/Projections/LeadingParticlesFinalState.hh"
9
10namespace Rivet {
11
12
13 /// ATLAS Z+jets in pp at 7 TeV
14 class ATLAS_2011_I945498 : public Analysis {
15 public:
16
17 /// Constructor
18 ATLAS_2011_I945498()
19 : Analysis("ATLAS_2011_I945498")
20 { }
21
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25
26 // Variable initialisation
27 _isZeeSample = false;
28 _isZmmSample = false;
29 for (size_t chn = 0; chn < 3; ++chn) {
30 book(weights_nj0[chn], "weights_nj0_" + to_str(chn));
31 book(weights_nj1[chn], "weights_nj1_" + to_str(chn));
32 book(weights_nj2[chn], "weights_nj2_" + to_str(chn));
33 book(weights_nj3[chn], "weights_nj3_" + to_str(chn));
34 book(weights_nj4[chn], "weights_nj4_" + to_str(chn));
35 }
36
37 // Set up projections
38 Cut cuts_mu = Cuts::abseta < 2.4 && Cuts::pT > 20*GeV;
39 DileptonFinder zfinder_mu(91.2*GeV, 0.1, cuts_mu && Cuts::abspid == PID::MUON, Cuts::massIn(66*GeV, 116*GeV));
40 declare(zfinder_mu, "DileptonFinder_mu");
41
42 Cut cuts_e = (Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47)) && Cuts::pT > 20*GeV;
43 DileptonFinder zfinder_el(91.2*GeV, 0.1, cuts_e && Cuts::abspid == PID::ELECTRON, Cuts::massIn(66*GeV, 116*GeV));
44 declare(zfinder_el, "DileptonFinder_el");
45
46 // For combined cross-sections (combined phase space + dressed level)
47 Cut cuts25_20 = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV;
48 DileptonFinder zfinder_comb_mu(91.2*GeV, 0.1, cuts25_20 && Cuts::abspid == PID::MUON, Cuts::massIn(66.0*GeV, 116.0*GeV));
49 declare(zfinder_comb_mu, "DileptonFinder_comb_mu");
50 DileptonFinder zfinder_comb_el(91.2*GeV, 0.1, cuts25_20 && Cuts::abspid == PID::ELECTRON, Cuts::massIn(66.0*GeV, 116.0*GeV));
51 declare(zfinder_comb_el, "DileptonFinder_comb_el");
52
53 // Define veto FS in order to prevent Z-decay products entering the jet algorithm
54 VetoedFinalState remfs;
55 remfs.addVetoOnThisFinalState(zfinder_el);
56 remfs.addVetoOnThisFinalState(zfinder_mu);
57 VetoedFinalState remfs_comb;
58 remfs_comb.addVetoOnThisFinalState(zfinder_comb_el);
59 remfs_comb.addVetoOnThisFinalState(zfinder_comb_mu);
60
61 FastJets jets(remfs, JetAlg::ANTIKT, 0.4);
62 jets.useInvisibles();
63 declare(jets, "jets");
64 FastJets jets_comb(remfs_comb, JetAlg::ANTIKT, 0.4);
65 jets_comb.useInvisibles();
66 declare(jets_comb, "jets_comb");
67
68 // 0=el, 1=mu, 2=comb
69 for (size_t chn = 0; chn < 3; ++chn) {
70 book(_h_njet_incl[chn] ,1, 1, chn+1);
71 book(_h_njet_ratio[chn] ,2, 1, chn+1);
72 book(_h_ptjet[chn] ,3, 1, chn+1);
73 book(_h_ptlead[chn] ,4, 1, chn+1);
74 book(_h_ptseclead[chn] ,5, 1, chn+1);
75 book(_h_yjet[chn] ,6, 1, chn+1);
76 book(_h_ylead[chn] ,7, 1, chn+1);
77 book(_h_yseclead[chn] ,8, 1, chn+1);
78 book(_h_mass[chn] ,9, 1, chn+1);
79 book(_h_deltay[chn] ,10, 1, chn+1);
80 book(_h_deltaphi[chn] ,11, 1, chn+1);
81 book(_h_deltaR[chn] ,12, 1, chn+1);
82 }
83 }
84
85
86 // Jet selection criteria universal for electron and muon channel
87 /// @todo Replace with a Cut passed to jetsByPt
88 Jets selectJets(const DileptonFinder* zf, const FastJets* allJets) {
89 const FourMomentum l1 = zf->constituents()[0].momentum();
90 const FourMomentum l2 = zf->constituents()[1].momentum();
91 Jets jets;
92 for (const Jet& jet : allJets->jetsByPt(Cuts::pT > 30*GeV)) {
93 const FourMomentum jmom = jet.momentum();
94 if (jmom.absrap() < 4.4 &&
95 deltaR(l1, jmom) > 0.5 && deltaR(l2, jmom) > 0.5) {
96 jets.push_back(jet);
97 }
98 }
99 return jets;
100 }
101
102
103 /// Perform the per-event analysis
104 void analyze(const Event& event) {
105 vector<const DileptonFinder*> zfs;
106 zfs.push_back(& (apply<DileptonFinder>(event, "DileptonFinder_el")));
107 zfs.push_back(& (apply<DileptonFinder>(event, "DileptonFinder_mu")));
108 zfs.push_back(& (apply<DileptonFinder>(event, "DileptonFinder_comb_el")));
109 zfs.push_back(& (apply<DileptonFinder>(event, "DileptonFinder_comb_mu")));
110
111 vector<const FastJets*> fjs;
112 fjs.push_back(& (apply<FastJets>(event, "jets")));
113 fjs.push_back(& (apply<FastJets>(event, "jets_comb")));
114
115 // Determine what kind of MC sample this is
116 const bool isZee = (zfs[0]->bosons().size() == 1) || (zfs[2]->bosons().size() == 1);
117 const bool isZmm = (zfs[1]->bosons().size() == 1) || (zfs[3]->bosons().size() == 1);
118 if (isZee) _isZeeSample = true;
119 if (isZmm) _isZmmSample = true;
120
121 // Require exactly one electronic or muonic Z-decay in the event
122 bool isZeemm = ( (zfs[0]->bosons().size() == 1 && zfs[1]->bosons().size() != 1) ||
123 (zfs[1]->bosons().size() == 1 && zfs[0]->bosons().size() != 1) );
124 bool isZcomb = ( (zfs[2]->bosons().size() == 1 && zfs[3]->bosons().size() != 1) ||
125 (zfs[3]->bosons().size() == 1 && zfs[2]->bosons().size() != 1) );
126 if (!isZeemm && !isZcomb) vetoEvent;
127
128 vector<int> zfIDs;
129 vector<int> fjIDs;
130 if (isZeemm) {
131 int chn = zfs[0]->bosons().size() == 1 ? 0 : 1;
132 zfIDs.push_back(chn);
133 fjIDs.push_back(0);
134 }
135 if (isZcomb) {
136 int chn = zfs[2]->bosons().size() == 1 ? 2 : 3;
137 zfIDs.push_back(chn);
138 fjIDs.push_back(1);
139 }
140
141 for (size_t izf = 0; izf < zfIDs.size(); ++izf) {
142 int zfID = zfIDs[izf];
143 int fjID = fjIDs[izf];
144
145 int chn = zfID;
146 if (zfID == 2 || zfID == 3) chn = 2;
147
148 Jets jets = selectJets(zfs[zfID], fjs[fjID]);
149
150 switch (jets.size()) {
151 case 0:
152 weights_nj0[chn]->fill();
153 break;
154 case 1:
155 weights_nj0[chn]->fill();
156 weights_nj1[chn]->fill();
157 break;
158 case 2:
159 weights_nj0[chn]->fill();
160 weights_nj1[chn]->fill();
161 weights_nj2[chn]->fill();
162 break;
163 case 3:
164 weights_nj0[chn]->fill();
165 weights_nj1[chn]->fill();
166 weights_nj2[chn]->fill();
167 weights_nj3[chn]->fill();
168 break;
169 default: // >= 4
170 weights_nj0[chn]->fill();
171 weights_nj1[chn]->fill();
172 weights_nj2[chn]->fill();
173 weights_nj3[chn]->fill();
174 weights_nj4[chn]->fill();
175 }
176
177 // Require at least one jet
178 if (jets.empty()) continue;
179
180 // Fill jet multiplicities
181 for (size_t ijet = 1; ijet <= jets.size(); ++ijet) {
182 _h_njet_incl[chn]->fill(ijet);
183 }
184
185 // Loop over selected jets, fill inclusive jet distributions
186 for (size_t ijet = 0; ijet < jets.size(); ++ijet) {
187 _h_ptjet[chn]->fill(jets[ijet].pT()/GeV);
188 _h_yjet [chn]->fill(fabs(jets[ijet].rapidity()));
189 }
190
191 // Leading jet histos
192 const double ptlead = jets[0].pT()/GeV;
193 const double yabslead = fabs(jets[0].rapidity());
194 _h_ptlead[chn]->fill(ptlead);
195 _h_ylead [chn]->fill(yabslead);
196
197 if (jets.size() >= 2) {
198 // Second jet histos
199 const double pt2ndlead = jets[1].pT()/GeV;
200 const double yabs2ndlead = fabs(jets[1].rapidity());
201 _h_ptseclead[chn] ->fill(pt2ndlead);
202 _h_yseclead [chn] ->fill(yabs2ndlead);
203
204 // Dijet histos
205 const double deltaphi = fabs(deltaPhi(jets[1], jets[0]));
206 const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity()) ;
207 const double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY));
208 const double mass = (jets[0].momentum() + jets[1].momentum()).mass();
209 _h_mass [chn] ->fill(mass/GeV);
210 _h_deltay [chn] ->fill(deltarap);
211 _h_deltaphi[chn] ->fill(deltaphi);
212 _h_deltaR [chn] ->fill(deltar);
213 }
214 }
215 }
216
217
218 /// @name Ratio calculator util functions
219 /// @{
220
221 /// Calculate the ratio, being careful about div-by-zero
222 double ratio(double a, double b) {
223 return (b != 0) ? a/b : 0;
224 }
225
226 /// Calculate the ratio error, being careful about div-by-zero
227 double ratio_err(double a, double b) {
228 return (b != 0) ? sqrt(a/sqr(b) + sqr(a)/(b*b*b)) : 0;
229 }
230
231 /// @}
232
233
234 void finalize() {
235 // Fill ratio histograms
236 for (size_t chn = 0; chn < 3; ++chn) {
237 _h_njet_ratio[chn]->bin(1).set(ratio(weights_nj1[chn]->val(), weights_nj0[chn]->val()),
238 ratio_err(weights_nj1[chn]->val(), weights_nj0[chn]->val()));
239 _h_njet_ratio[chn]->bin(2).set(ratio(weights_nj2[chn]->val(), weights_nj1[chn]->val()),
240 ratio_err(weights_nj2[chn]->val(), weights_nj1[chn]->val()));
241 _h_njet_ratio[chn]->bin(3).set(ratio(weights_nj3[chn]->val(), weights_nj2[chn]->val()),
242 ratio_err(weights_nj3[chn]->val(), weights_nj2[chn]->val()));
243 _h_njet_ratio[chn]->bin(4).set(ratio(weights_nj4[chn]->val(), weights_nj3[chn]->val()),
244 ratio_err(weights_nj4[chn]->val(), weights_nj3[chn]->val()));
245 }
246
247 // Scale other histos
248 for (size_t chn = 0; chn < 3; ++chn) {
249 // For ee and mumu channels: normalize to Njet inclusive cross-section
250 double xs = crossSectionPerEvent()/picobarn;
251 if (chn != 2 && weights_nj0[chn]->val() != 0.) xs = 1.0 / weights_nj0[chn]->val();
252 // For inclusive MC sample(ee/mmu channels together) we want the single-lepton-flavor xsec
253 if (_isZeeSample && _isZmmSample) xs *= 0.5;
254
255 // Special case histogram: always not normalized
256 scale(_h_njet_incl[chn], (chn < 2) ? crossSectionPerEvent()/picobarn : xs);
257
258 scale(_h_ptjet[chn] , xs);
259 scale(_h_ptlead[chn] , xs);
260 scale(_h_ptseclead[chn], xs);
261 scale(_h_yjet[chn] , xs);
262 scale(_h_ylead[chn] , xs);
263 scale(_h_yseclead[chn] , xs);
264 scale(_h_deltaphi[chn] , xs);
265 scale(_h_deltay[chn] , xs);
266 scale(_h_deltaR[chn] , xs);
267 scale(_h_mass[chn] , xs);
268 }
269
270 }
271
272 /// @}
273
274
275 private:
276
277 bool _isZeeSample;
278 bool _isZmmSample;
279
280 CounterPtr weights_nj0[3];
281 CounterPtr weights_nj1[3];
282 CounterPtr weights_nj2[3];
283 CounterPtr weights_nj3[3];
284 CounterPtr weights_nj4[3];
285
286 Estimate1DPtr _h_njet_ratio[3];
287 Histo1DPtr _h_njet_incl[3];
288 Histo1DPtr _h_ptjet[3];
289 Histo1DPtr _h_ptlead[3];
290 Histo1DPtr _h_ptseclead[3];
291 Histo1DPtr _h_yjet[3];
292 Histo1DPtr _h_ylead[3];
293 Histo1DPtr _h_yseclead[3];
294 Histo1DPtr _h_deltaphi[3];
295 Histo1DPtr _h_deltay[3];
296 Histo1DPtr _h_deltaR[3];
297 Histo1DPtr _h_mass[3];
298
299 };
300
301
302 RIVET_DECLARE_PLUGIN(ATLAS_2011_I945498);
303
304
305}
|