Rivet analyses referenceATLAS_2013_I1230812$Z$ + jets in $pp$ at 7 TeVExperiment: ATLAS (LHC) Inspire ID: 1230812 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
Measurements of the production of jets of particles in association with a $Z$ boson in $pp$ collisions at $\sqrt{s}$ = 7 TeV are presented, using data corresponding to an integrated luminosity of 4.6/fb collected by the ATLAS experiment at the Large Hadron Collider. Inclusive and differential jet cross sections in $Z$ events, with Z decaying into electron or muon pairs, are measured for jets with transverse momentum $p_T > 30$ GeV and rapidity $|y| < 4.4$. This Rivet module implements the event selection for the weighted combination of both decay channels to produce the average cross section for a single lepton flavour, and uses the data from that combination (as in the paper plots). In the default mode this is how it will run, assuming mixed electronic and muonic events. If LMODE is set to EL (MU), the plots for the individual decay channels (with a slightly different fiducial phase space) are made. Source code: ATLAS_2013_I1230812.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/DileptonFinder.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/VetoedFinalState.hh"
6#include "Rivet/Projections/PromptFinalState.hh"
7
8namespace Rivet {
9
10
11 /// Z + jets in pp at 7 TeV (combined channel / base class)
12 ///
13 /// @note This base class contains a "mode" variable for combined, e, and mu channel derived classes
14 class ATLAS_2013_I1230812 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2013_I1230812);
19
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // Get options from the new option system
25 _mode = 0;
26 if ( getOption("LMODE") == "EL" ) _mode = 1;
27 if ( getOption("LMODE") == "MU" ) _mode = 2;
28
29 // Determine the e/mu decay channels used (NB Prompt leptons only).
30 /// @todo Note that Zs are accepted with any rapidity: the cuts are on the e/mu: is this correct?
31 Cut pt20 = Cuts::pT >= 20*GeV;
32 Cut eta_e = _mode? Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47) : Cuts::abseta < 2.5;
33 Cut eta_m = _mode? Cuts::abseta < 2.4 : Cuts::abseta < 2.5;
34 DileptonFinder zfinder_el(91.2*GeV, 0.1, eta_e && pt20 && Cuts::abspid == PID::ELECTRON, Cuts::massIn(66*GeV, 116*GeV));
35 DileptonFinder zfinder_mu(91.2*GeV, 0.1, eta_m && pt20 && Cuts::abspid == PID::MUON, Cuts::massIn(66*GeV, 116*GeV));
36 declare(zfinder_el, "zfinder_el");
37 declare(zfinder_mu, "zfinder_mu");
38
39 // Define veto FS in order to prevent Z-decay products entering the jet algorithm
40 VetoedFinalState had_fs;
41 had_fs.addVetoOnThisFinalState(getProjection<DileptonFinder>("zfinder_el"));
42 had_fs.addVetoOnThisFinalState(getProjection<DileptonFinder>("zfinder_mu"));
43 FastJets jets(had_fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL);
44 declare(jets, "jets");
45
46 book(_h_njet_incl , 1, 1, _mode+1);
47 book(_h_njet_incl_ratio , 2, 1, _mode+1);
48 book(_h_njet_excl , 3, 1, _mode+1);
49 book(_h_njet_excl_ratio , 4, 1, _mode+1);
50 book(_h_njet_excl_pt150 , 5, 1, _mode+1);
51 book(_h_njet_excl_pt150_ratio , 6, 1, _mode+1);
52 book(_h_njet_excl_vbf , 7, 1, _mode+1);
53 book(_h_njet_excl_vbf_ratio , 8, 1, _mode+1);
54 book(_h_ptlead , 9, 1, _mode+1);
55 book(_h_ptseclead , 10, 1, _mode+1);
56 book(_h_ptthirdlead , 11, 1, _mode+1);
57 book(_h_ptfourthlead , 12, 1, _mode+1);
58 book(_h_ptlead_excl , 13, 1, _mode+1);
59 book(_h_pt_ratio , 14, 1, _mode+1);
60 book(_h_pt_z , 15, 1, _mode+1);
61 book(_h_pt_z_excl , 16, 1, _mode+1);
62 book(_h_ylead , 17, 1, _mode+1);
63 book(_h_yseclead , 18, 1, _mode+1);
64 book(_h_ythirdlead , 19, 1, _mode+1);
65 book(_h_yfourthlead , 20, 1, _mode+1);
66 book(_h_deltay , 21, 1, _mode+1);
67 book(_h_mass , 22, 1, _mode+1);
68 book(_h_deltaphi , 23, 1, _mode+1);
69 book(_h_deltaR , 24, 1, _mode+1);
70 book(_h_ptthirdlead_vbf , 25, 1, _mode+1);
71 book(_h_ythirdlead_vbf , 26, 1, _mode+1);
72 book(_h_ht , 27, 1, _mode+1);
73 book(_h_st , 28, 1, _mode+1);
74 }
75
76
77 /// Perform the per-event analysis
78 void analyze(const Event& event) {
79
80 FourMomentum z, lp, lm;
81 const DileptonFinder& zfinder_el = apply<DileptonFinder>(event, "zfinder_el");
82 const DileptonFinder& zfinder_mu = apply<DileptonFinder>(event, "zfinder_mu");
83
84 bool e_ok = zfinder_el.constituents().size() == 2 && zfinder_mu.constituents().size() == 0;
85 bool m_ok = zfinder_el.constituents().size() == 0 && zfinder_mu.constituents().size() == 2;
86
87 if (_mode == 0 && !e_ok && !m_ok) vetoEvent;
88 if (_mode == 1 && !e_ok) vetoEvent;
89 if (_mode == 2 && !m_ok) vetoEvent;
90
91 if (zfinder_el.constituents().size() == 2) {
92 z = zfinder_el.boson().momentum();
93 lp = zfinder_el.constituents()[0].momentum();
94 lm = zfinder_el.constituents()[1].momentum();
95 }
96 else if (zfinder_mu.constituents().size() == 2) {
97 z = zfinder_mu.boson().momentum();
98 lp = zfinder_mu.constituents()[0].momentum();
99 lm = zfinder_mu.constituents()[1].momentum();
100 }
101 else vetoEvent;
102
103 if (deltaR(lp, lm) < 0.2) vetoEvent;
104
105 Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4);
106 idiscard(jets, deltaRLess(lp, 0.5));
107 idiscard(jets, deltaRLess(lm, 0.5));
108
109 // Fill jet multiplicities
110 for (size_t ijet = 0; ijet <= jets.size(); ++ijet) {
111 _h_njet_incl->fill(ijet);
112 }
113 _h_njet_excl->fill(jets.size());
114
115 // Require at least one jet
116 if (jets.size() >= 1) {
117 // Leading jet histos
118 const double ptlead = jets[0].pT()/GeV;
119 const double yabslead = fabs(jets[0].rapidity());
120 const double ptz = z.pT()/GeV;
121 _h_ptlead->fill(ptlead);
122 _h_ylead ->fill(yabslead);
123 _h_pt_z ->fill(ptz);
124 // Fill jet multiplicities
125 if (ptlead > 150) _h_njet_excl_pt150->fill(jets.size());
126
127 // Loop over selected jets, fill inclusive distributions
128 double st = 0;
129 double ht = lp.pT()/GeV + lm.pT()/GeV;
130 for (size_t ijet = 0; ijet < jets.size(); ++ijet) {
131 ht += jets[ijet].pT()/GeV;
132 st += jets[ijet].pT()/GeV;
133 }
134 _h_ht->fill(ht);
135 _h_st->fill(st);
136
137 // Require exactly one jet
138 if (jets.size() == 1) {
139 _h_ptlead_excl->fill(ptlead);
140 _h_pt_z_excl ->fill(ptz);
141 }
142 }
143
144
145 // Require at least two jets
146 if (jets.size() >= 2) {
147 // Second jet histos
148 const double ptlead = jets[0].pT()/GeV;
149 const double pt2ndlead = jets[1].pT()/GeV;
150 const double ptratio = pt2ndlead/ptlead;
151 const double yabs2ndlead = fabs(jets[1].rapidity());
152 _h_ptseclead->fill(pt2ndlead);
153 _h_yseclead->fill( yabs2ndlead);
154 _h_pt_ratio->fill( ptratio);
155
156 // Dijet histos
157 const double deltaphi = fabs(deltaPhi(jets[1], jets[0]));
158 const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity()) ;
159 const double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY));
160 const double mass = (jets[0].momentum() + jets[1].momentum()).mass()/GeV;
161 _h_mass->fill( mass);
162 _h_deltay->fill( deltarap);
163 _h_deltaphi->fill(deltaphi);
164 _h_deltaR->fill( deltar);
165
166 if (mass > 350 && deltarap > 3) _h_njet_excl_vbf->fill(jets.size());
167 }
168
169 // Require at least three jets
170 if (jets.size() >= 3) {
171 // Third jet histos
172 const double pt3rdlead = jets[2].pT()/GeV;
173 const double yabs3rdlead = fabs(jets[2].rapidity());
174 _h_ptthirdlead->fill(pt3rdlead);
175 _h_ythirdlead->fill( yabs3rdlead);
176
177 //Histos after VBF preselection
178 const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity()) ;
179 const double mass = (jets[0].momentum() + jets[1].momentum()).mass();
180 if (mass > 350 && deltarap > 3) {
181 _h_ptthirdlead_vbf->fill(pt3rdlead);
182 _h_ythirdlead_vbf->fill( yabs3rdlead);
183 }
184 }
185
186 // Require at least four jets
187 if (jets.size() >= 4) {
188 // Fourth jet histos
189 const double pt4thlead = jets[3].pT()/GeV;
190 const double yabs4thlead = fabs(jets[3].rapidity());
191 _h_ptfourthlead->fill(pt4thlead);
192 _h_yfourthlead->fill( yabs4thlead);
193 }
194 }
195
196 /// @name Ratio calculator util functions
197 /// @{
198
199 /// Calculate the efficiency error, being careful about div-by-zero
200 double err_incl(const YODA::Dbn1D& M, const YODA::Dbn1D& N, bool hasWeights) {
201 double r = safediv(M.sumW(), N.sumW());
202 if (hasWeights) { // use F. James's approximation for weighted events
203 return sqrt( safediv((1 - 2 * r) * M.sumW2() + r * r * N.sumW2(), N.sumW() * N.sumW()) );
204 }
205 return sqrt( safediv(r * (1 - r), N.sumW()) );
206 }
207
208 /// Calculate the ratio error, being careful about div-by-zero
209 double err_excl(const YODA::Dbn1D& A, const YODA::Dbn1D& B) {
210 double r = safediv(A.sumW(), B.sumW());
211 double dAsquared = safediv(A.sumW2(), A.sumW() * A.sumW()); // squared relative error of A
212 double dBsquared = safediv(B.sumW2(), B.sumW() * B.sumW()); // squared relative error of B
213 return r * sqrt(dAsquared + dBsquared);
214 }
215
216 /// @}
217
218
219 void finalize() {
220 const bool hasWeights = _h_njet_incl->effNumEntries() != _h_njet_incl->numEntries();
221 for (size_t i = 0; i < 6; ++i) {
222 _h_njet_incl_ratio->bin(i+1).set(safediv(_h_njet_incl->bin(i + 2).sumW(), _h_njet_incl->bin(i+1).sumW()),
223 err_incl(_h_njet_incl->bin(i + 2).raw(), _h_njet_incl->bin(i+1).raw(), hasWeights));
224 _h_njet_excl_ratio->bin(i+1).set(safediv(_h_njet_excl->bin(i + 2).sumW(), _h_njet_excl->bin(i+1).sumW()),
225 err_excl(_h_njet_excl->bin(i + 2).raw(), _h_njet_excl->bin(i+1).raw()));
226 if (i >= 1) {
227 _h_njet_excl_pt150_ratio->bin(i).set(safediv(_h_njet_excl_pt150->bin(i+1).sumW(), _h_njet_excl_pt150->bin(i).sumW()),
228 err_excl(_h_njet_excl_pt150->bin(i+1).raw(), _h_njet_excl_pt150->bin(i).raw()));
229 if (i >= 2) {
230 _h_njet_excl_vbf_ratio->bin(i - 1).set(safediv(_h_njet_excl_vbf->bin(i+1).sumW(), _h_njet_excl_vbf->bin(i).sumW()),
231 err_excl(_h_njet_excl_vbf->bin(i+1).raw(), _h_njet_excl_vbf->bin(i).raw()));
232 }
233 }
234 }
235
236 const double sf = _mode? 1.0 : 0.5;
237 const double xs = sf * crossSectionPerEvent()/picobarn;
238
239 scale(_h_njet_incl, xs); scale(_h_njet_excl, xs); scale(_h_njet_excl_pt150, xs);
240 scale(_h_njet_excl_vbf, xs); scale(_h_ptlead, xs); scale(_h_ptseclead, xs);
241 scale(_h_ptthirdlead, xs); scale(_h_ptfourthlead, xs); scale(_h_ptlead_excl, xs);
242 scale(_h_pt_ratio, xs); scale(_h_pt_z, xs); scale(_h_pt_z_excl, xs);
243 scale(_h_ylead, xs); scale(_h_yseclead, xs); scale(_h_ythirdlead, xs);
244 scale(_h_yfourthlead, xs); scale(_h_deltay, xs); scale(_h_mass, xs);
245 scale(_h_deltaphi, xs); scale(_h_deltaR, xs); scale(_h_ptthirdlead_vbf, xs);
246 scale(_h_ythirdlead_vbf, xs); scale(_h_ht, xs); scale(_h_st, xs);
247 }
248
249 /// @}
250
251
252 protected:
253
254 size_t _mode;
255
256
257 private:
258
259 Estimate1DPtr _h_njet_incl_ratio;
260 Estimate1DPtr _h_njet_excl_ratio;
261 Estimate1DPtr _h_njet_excl_pt150_ratio;
262 Estimate1DPtr _h_njet_excl_vbf_ratio;
263 Histo1DPtr _h_njet_incl;
264 Histo1DPtr _h_njet_excl;
265 Histo1DPtr _h_njet_excl_pt150;
266 Histo1DPtr _h_njet_excl_vbf;
267 Histo1DPtr _h_ptlead;
268 Histo1DPtr _h_ptseclead;
269 Histo1DPtr _h_ptthirdlead;
270 Histo1DPtr _h_ptfourthlead;
271 Histo1DPtr _h_ptlead_excl;
272 Histo1DPtr _h_pt_ratio;
273 Histo1DPtr _h_pt_z;
274 Histo1DPtr _h_pt_z_excl;
275 Histo1DPtr _h_ylead;
276 Histo1DPtr _h_yseclead;
277 Histo1DPtr _h_ythirdlead;
278 Histo1DPtr _h_yfourthlead;
279 Histo1DPtr _h_deltay;
280 Histo1DPtr _h_mass;
281 Histo1DPtr _h_deltaphi;
282 Histo1DPtr _h_deltaR;
283 Histo1DPtr _h_ptthirdlead_vbf;
284 Histo1DPtr _h_ythirdlead_vbf;
285 Histo1DPtr _h_ht;
286 Histo1DPtr _h_st;
287 };
288
289
290 RIVET_DECLARE_PLUGIN(ATLAS_2013_I1230812);
291
292}
|
