rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2013_I1230812

$Z$ + jets in $pp$ at 7 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1230812
Status: VALIDATED
Authors:
  • Katharina Bierwagen
  • Frank Siegert
References: Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • Z+jets, electronic Z-decays (data are a weighted combination of electron/muon).

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}