rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

LHCB_2016_I1454404

Measurement of forward W and Z boson production in association with jets at LHCb
Experiment: LHCB (LHC)
Inspire ID: 1454404
Status: VALIDATED
Authors:
  • Abbie Jane Chadwick
  • Stephen Farry
References: Beams: p+ p+
Beam energies: (4000.0, 4000.0) GeV
Run details:
  • p+ p+ -> W/Z + jet + X

Measurements are made of forward W and Z production in association with jets in the forward region. Muons and jets are required to have transverse momentum in excess of 20 Gev, and to have a pseudorapidity between 2.0 and 4.5 for muons, and between 2.2 and 4.2 for jets. A single muon is required in the case of $W$ production, and two opposite sign muons with a combined invariant mass of between 60 and 120 GeV are required for $W$ production. The leptons and jets are required to be separated by a radius of 0.5 in ($\eta$, $\phi$) space.

Source code: LHCB_2016_I1454404.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/DileptonFinder.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Measurement of forward W and Z boson production with jets in pp collisions at 8 TeV
 13  class LHCB_2016_I1454404 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2016_I1454404);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      _mode = 0;
 27      string mode = getOption("MODE");
 28      if (mode == "ALL" ) _mode = 0;
 29      else if (mode == "WpJET") _mode = 1;
 30      else if (mode == "WmJET") _mode = 2;
 31      else if (mode == "ZJET") _mode = 3;
 32      else if (mode == "WJET") _mode = 4;
 33
 34      const Cut muSel = Cuts::eta >= 2.0 && Cuts::eta <= 4.5 && Cuts::pT > 20*GeV;
 35
 36      // Z boson
 37      DileptonFinder zfinder(91.2*GeV, 0.1, muSel && Cuts::abspid == PID::MUON, Cuts::massIn(60*GeV, 120*GeV));
 38      declare(zfinder, "DileptonFinder");
 39
 40      // W boson
 41      declare("MET", MissingMomentum());
 42      LeptonFinder mf(muSel && Cuts::abspid == PID::MUON, 0.0);
 43      declare(mf, "Muons");
 44
 45      // Jet Z
 46      FastJets jetproZ(zfinder.remainingFinalState(), JetAlg::ANTIKT, 0.5); ///< @todo No jet acceptance?
 47      declare(jetproZ, "JetsZ");
 48
 49      // Jet W
 50      VetoedFinalState nomufs; ///< @todo No jet acceptance?
 51      nomufs.vetoFinalState(mf);
 52      FastJets jetproW(nomufs, JetAlg::ANTIKT, 0.5);
 53      declare(jetproW, "JetsW");
 54
 55
 56      // Book histograms
 57      if (_mode == 0 || _mode == 1 || _mode == 4) {
 58        book(_h_wpj, 1, 1, 1);
 59        book(_h_eta_wpj, 4, 1, 1);
 60        book(_h_etaj_wpj, 5, 1, 1);
 61        book(_h_ptj_wpj, 6, 1, 1);
 62      }
 63      if (_mode == 0 || _mode == 2 || _mode == 4) {
 64        book(_h_wmj, 1, 1, 2);
 65        book(_h_eta_wmj, 4, 1, 2);
 66        book(_h_etaj_wmj, 5, 1, 2);
 67        book(_h_ptj_wmj, 6, 1, 2);
 68      }
 69      if (_mode == 0 || _mode == 3) {
 70        book(_h_zj, 1, 1, 3);
 71        book(_h_yz_zj, 7, 1, 1);
 72        book(_h_etaj_zj, 8, 1, 1);
 73        book(_h_ptj_zj, 9, 1, 1);
 74        book(_h_dphi_zj, 10, 1, 1);
 75      }
 76      if (_mode == 0 ){
 77        book(_h_rwz, 2,1,1);
 78        book(_h_rwpz, 2,1,2);
 79        book(_h_rwmz, 2,1,3);
 80      }
 81      if (_mode == 0 || _mode == 4){
 82        book(_h_rwpm, 2,1,4);
 83        book(_h_aw, 3,1,1);
 84        // this is a temporary histogram to construct rwz later
 85        book(_h_wj, "_temp_wj", refData(1,1,1));
 86      }
 87    }
 88
 89
 90    /// Perform the per-event analysis
 91    void analyze(const Event& event) {
 92      const Cut jetSel = Cuts::eta >= 2.2 && Cuts::eta <= 4.2 && Cuts::pT > 20*GeV;
 93
 94      // Z
 95      if (_mode == 0 || _mode == 3) {
 96
 97        const DileptonFinder& zfinder = apply<DileptonFinder>(event, "DileptonFinder");
 98        if (zfinder.bosons().size() ==1){
 99          const Particles muon = zfinder.leptons();
100          const Particles Z = zfinder.bosons();
101          const FourMomentum Zmom = Z[0].momentum();
102          const Jets jetsZ = apply<FastJets>(event, "JetsZ").jetsByPt(jetSel);
103          const Jets cleanedJetsZ = discard(jetsZ, [&](const Jet& j) {return any(muon, deltaRLess(j, 0.5)); });
104
105          if (cleanedJetsZ.size() > 0 && cleanedJetsZ.at(0).pT() > 20*GeV) {
106            const double yZ = Zmom.rap(); //histogram 7
107            const double etaj = cleanedJetsZ[0].eta(); //histogram 8
108            const double ptj  = cleanedJetsZ[0].pT()/GeV; //histogram 9
109            double dphi_tmp = abs(Zmom.phi() - cleanedJetsZ[0].phi());
110            const double dphi = dphi_tmp < Rivet::pi ? dphi_tmp : Rivet::twopi - dphi_tmp;
111            _h_zj->fill(sqrtS()/GeV);
112            _h_dphi_zj->fill(dphi);
113            _h_yz_zj->fill(yZ); // boson rapidity vs diff cross section
114            _h_etaj_zj->fill(etaj); // jet pseudorapidity vs diff cross section
115            _h_ptj_zj->fill(ptj); //jet transverse momentum vs diff cross section
116          }
117        }
118      }
119
120      // W
121      if (_mode == 0 || _mode == 1 || _mode == 2  || _mode == 4) {
122
123        // Identify the closest-matching l+MET to m == mW
124        const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
125        const Particles& mus = apply<LeptonFinder>(event, "Muons").particles();
126        const int ifound = closestMatchIndex(mus, pmiss, Kin::mass, 80.4*GeV, 0*GeV, 500*GeV);
127        if (ifound >= 0) {
128          const Particle& mu = mus[ifound];
129
130          const Jets jetsW = apply<FastJets>(event, "JetsW").jetsByPt(jetSel);
131          const Jets cleanedJetsW = discard(jetsW, deltaRLess(mu, 0.5));
132
133          if (cleanedJetsW.size() > 0 && cleanedJetsW.at(0).pT() > 20*GeV) {
134            const double etaj = cleanedJetsW[0].eta(); //histogram 5
135            const double etamu = mu.eta(); //histogram 4
136            if ( (_mode == 0 || _mode == 1 || _mode == 4) && mu.charge() > 0) {
137              //fill with W related analysis
138              if (_mode != 1 ) _h_wj->fill(sqrtS()/GeV); // don't need this for single charge case
139              _h_wpj->fill(sqrtS()/GeV);
140              _h_eta_wpj->fill(etamu); // W+ Jet  muon pseudorapidity vs diff cross section
141              _h_etaj_wpj->fill(etaj); // W+ Jet jet pseudorapidity vs diff cross section
142              _h_ptj_wpj->fill(cleanedJetsW[0].pT()/GeV); // W+ Jet jet transverse momentum vs diff cross section
143            }
144            else if ( (_mode == 0 || _mode == 2 || _mode == 4 ) && mu.charge() < 0) {
145              //fill with W related analysis
146              if (_mode != 2) _h_wj->fill(sqrtS()/GeV);
147              _h_wmj->fill(sqrtS()/GeV); // don't need this for single charge case
148              _h_eta_wmj->fill(etamu); // W- Jet  muon pseudorapidity vs diff cross section
149              _h_etaj_wmj->fill(etaj); // W- Jet jet pseudorapidity vs diff cross section
150              _h_ptj_wmj->fill(cleanedJetsW[0].pT()/GeV); // W+ Jet jet transverse momentum vs diff cross section
151            }
152          }
153        }
154
155      }
156    }
157
158
159    /// Normalise histograms etc., after the run
160    void finalize() {
161      const double scalefactor = crossSection()/picobarn/sumOfWeights();
162
163      if (_mode == 0 || _mode == 1 || _mode == 4) {
164        scale({_h_wpj, _h_eta_wpj, _h_etaj_wpj, _h_ptj_wpj}, scalefactor);
165      }
166      if (_mode == 0 || _mode == 2 || _mode == 4) {
167        scale({_h_wmj, _h_eta_wmj, _h_etaj_wmj, _h_ptj_wmj}, scalefactor);
168      }
169      if (_mode == 0 || _mode == 3) {
170        scale({_h_zj, _h_yz_zj, _h_etaj_zj, _h_ptj_zj, _h_dphi_zj}, scalefactor);
171      }
172      if (_mode == 0 ) {
173        scale(_h_wj, scalefactor); // need to scale this for consistency
174        divide(_h_wpj, _h_zj, _h_rwpz);
175        divide(_h_wmj, _h_zj, _h_rwmz);
176        divide(_h_wj, _h_zj, _h_rwz);
177      }
178      if (_mode == 0 || _mode == 4) {
179        divide(_h_wpj, _h_wmj, _h_rwpm);
180        asymm(_h_wpj, _h_wmj, _h_aw);
181      }
182    }
183
184
185  protected:
186
187    size_t _mode;
188
189    /// @name Histograms
190    /// @{
191    Histo1DPtr _h_wpj, _h_wmj, _h_wj, _h_zj;
192    Estimate1DPtr _h_rwz, _h_rwpz, _h_rwmz, _h_rwpm, _h_aw;
193    Histo1DPtr _h_eta_wpj, _h_eta_wmj, _h_etaj_wpj, _h_etaj_wmj, _h_ptj_wpj, _h_ptj_wmj;
194    Histo1DPtr _h_yz_zj, _h_etaj_zj, _h_ptj_zj, _h_dphi_zj;
195    /// @}
196
197  };
198
199
200  RIVET_DECLARE_PLUGIN(LHCB_2016_I1454404);
201
202}