rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2014_I1304688

Measurement of jet multiplicity and transverse momentum spectra in top events using full 7 TeV ATLAS dataset
Experiment: ATLAS (LHC)
Inspire ID: 1304688
Status: VALIDATED
Authors:
  • W. H. Bell
  • A. Grohsjean
References: Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • ttbar events with at least one lepton in the ttbar final state at 7 TeV, i.e. both semileptonic and dileptonic decays should be enabled. The tau decay channels also count as leptonic.

Measurement of the differential t$\bar{\mathrm t}$ production cross-section in 7 TeV proton-proton collisions in the single-lepton channel from ATLAS. The data comprise the full 2011 data sample corresponding to an integrated luminosity of $4.6\,\mathrm{fb}^{-1}$. The differential cross-sections are measured as a function of the jet multiplicity for up to eight jets using jet transverse momentum thresholds of 25, 40, 60, and 80 GeV, and as a function of jet transverse momentum up to the fifth leading jet. The results after background subtraction are corrected for all detector effects, within a kinematic range closely matched to the experimental acceptance.

Source code: ATLAS_2014_I1304688.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/VetoedFinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/DressedLeptons.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include <bitset>
  9
 10namespace Rivet {
 11
 12
 13  /// @brief ATLAS 7 TeV jets in ttbar events analysis
 14  ///
 15  /// @author W. H. Bell <W.Bell@cern.ch>
 16  /// @author A. Grohsjean <alexander.grohsjean@desy.de>
 17  class ATLAS_2014_I1304688 : public Analysis {
 18  public:
 19
 20    ATLAS_2014_I1304688():
 21      Analysis("ATLAS_2014_I1304688"),
 22      _jet_ntag(0),
 23      _met_et(0.),
 24      _met_phi(0.),
 25      _hMap(),
 26      //_chanLimit(3),
 27      _histLimit(6)
 28    {   }
 29
 30
 31    void init() {
 32      // Eta ranges
 33      /// @todo 1 MeV? Really?
 34      Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV;
 35      Cut eta_lep = Cuts::abseta < 2.5;
 36
 37      // All final state particles
 38      FinalState fs(eta_full);
 39
 40      // Get photons to dress leptons
 41      IdentifiedFinalState photons(fs);
 42      photons.acceptIdPair(PID::PHOTON);
 43
 44      // Projection to find the electrons
 45      IdentifiedFinalState el_id(fs);
 46      el_id.acceptIdPair(PID::ELECTRON);
 47      PromptFinalState electrons(el_id);
 48      electrons.acceptTauDecays(true);
 49      declare(electrons, "electrons");
 50      DressedLeptons dressedelectrons(photons, electrons, 0.1, eta_lep && Cuts::pT > 25*GeV, true);
 51      declare(dressedelectrons, "dressedelectrons");
 52      DressedLeptons vetodressedelectrons(photons, electrons, 0.1, eta_lep && Cuts::pT >= 15*GeV, true);
 53      declare(vetodressedelectrons, "vetodressedelectrons");
 54      DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full, true);
 55      declare(ewdressedelectrons, "ewdressedelectrons");
 56
 57      // Projection to find the muons
 58      IdentifiedFinalState mu_id(fs);
 59      mu_id.acceptIdPair(PID::MUON);
 60      PromptFinalState muons(mu_id);
 61      muons.acceptTauDecays(true);
 62      declare(muons, "muons");
 63      vector<pair<double, double> > eta_muon;
 64      DressedLeptons dressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT >= 25*GeV, true);
 65      declare(dressedmuons, "dressedmuons");
 66      DressedLeptons vetodressedmuons(photons, muons, 0.1, eta_lep && Cuts::pT >= 15*GeV, true);
 67      declare(vetodressedmuons, "vetodressedmuons");
 68      DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full, true);
 69      declare(ewdressedmuons, "ewdressedmuons");
 70
 71      // Projection to find neutrinos and produce MET
 72      IdentifiedFinalState nu_id;
 73      nu_id.acceptNeutrinos();
 74      PromptFinalState neutrinos(nu_id);
 75      neutrinos.acceptTauDecays(true);
 76      declare(neutrinos, "neutrinos");
 77
 78      // Jet clustering.
 79      VetoedFinalState vfs;
 80      vfs.addVetoOnThisFinalState(ewdressedelectrons);
 81      vfs.addVetoOnThisFinalState(ewdressedmuons);
 82      vfs.addVetoOnThisFinalState(neutrinos);
 83      FastJets jets(vfs, FastJets::ANTIKT, 0.4);
 84      jets.useInvisibles();
 85      declare(jets, "jets");
 86
 87      // Book histograms
 88      for (unsigned int ihist = 0; ihist < _histLimit ; ihist++) {
 89        const unsigned int threshLimit = _thresholdLimit(ihist);
 90        for (unsigned int ithres = 0; ithres < threshLimit; ithres++) {
 91          _histogram(ihist, ithres); // Create all histograms
 92        }
 93      }
 94    }
 95
 96
 97    void analyze(const Event& event) {
 98
 99      // Get the selected objects, using the projections.
100      _dressedelectrons = sortByPt(apply<DressedLeptons>(event, "dressedelectrons").dressedLeptons());
101      _vetodressedelectrons = apply<DressedLeptons>(event, "vetodressedelectrons").dressedLeptons();
102
103      _dressedmuons = sortByPt(apply<DressedLeptons>(event, "dressedmuons").dressedLeptons());
104      _vetodressedmuons = apply<DressedLeptons>(event, "vetodressedmuons").dressedLeptons();
105
106      _neutrinos = apply<PromptFinalState>(event, "neutrinos").particlesByPt();
107
108      _jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
109
110
111      // Calculate the missing ET, using the prompt neutrinos only (really?)
112      /// @todo Why not use MissingMomentum?
113      FourMomentum pmet;
114      for (const Particle& p : _neutrinos) pmet += p.momentum();
115      _met_et = pmet.pT();
116      _met_phi = pmet.phi();
117
118      // Check overlap of jets/leptons.
119      unsigned int i,j;
120      _jet_ntag = 0;
121      _overlap = false;
122      for (i = 0; i < _jets.size(); i++) {
123        const Jet& jet = _jets[i];
124        // If dR(el,jet) < 0.4 skip the event
125        for (const DressedLepton& el : _dressedelectrons) {
126          if (deltaR(jet, el) < 0.4) _overlap = true;
127        }
128        // If dR(mu,jet) < 0.4 skip the event
129        for (const DressedLepton& mu : _dressedmuons) {
130          if (deltaR(jet, mu) < 0.4) _overlap = true;
131        }
132        // If dR(jet,jet) < 0.5 skip the event
133        for (j = 0; j < _jets.size(); j++) {
134          const Jet& jet2 = _jets[j];
135          if (i == j) continue; // skip the diagonal
136          if (deltaR(jet, jet2) < 0.5) _overlap = true;
137        }
138        // Count the number of b-tags
139        if (!jet.bTags().empty()) _jet_ntag += 1;
140      }
141
142      // Evaluate basic event selection
143      unsigned int ejets_bits = 0, mujets_bits = 0;
144      bool pass_ejets = _ejets(ejets_bits);
145      bool pass_mujets = _mujets(mujets_bits);
146
147      // Remove events with object overlap
148      if (_overlap) vetoEvent;
149      // basic event selection requirements
150      if (!pass_ejets && !pass_mujets) vetoEvent;
151
152      // Check if the additional pT threshold requirements are passed
153      bool pass_jetPt = _additionalJetCuts();
154
155      // Count the jet multiplicity for 25, 40, 60 and 80GeV
156      unsigned int ithres, jet_n[4];
157      for (ithres = 0; ithres < 4; ithres++) {
158        jet_n[ithres] = _countJets(ithres);
159      }
160
161      // Fill histograms
162      for (unsigned int ihist = 0; ihist < 6; ihist++) {
163        if (ihist > 0 && !pass_jetPt) continue; // additional pT threshold cuts for pT plots
164        unsigned int threshLimit = _thresholdLimit(ihist);
165        for (ithres = 0; ithres < threshLimit; ithres++) {
166          if (jet_n[ithres] < 3) continue; // 3 or more jets for ljets
167          // Fill
168          if (ihist == 0) _histogram(ihist, ithres)->fill(jet_n[ithres]-2); // njets
169          else if (ihist == 1) _histogram(ihist, ithres)->fill(_jets[0].pT()); // leading jet pT
170          else if (ihist == 2) _histogram(ihist, ithres)->fill(_jets[1].pT()); // 2nd jet pT
171          else if (ihist == 3 && jet_n[ithres] >= 3) _histogram(ihist, ithres)->fill(_jets[2].pT()); // 3rd jet pT
172          else if (ihist == 4 && jet_n[ithres] >= 4) _histogram(ihist, ithres)->fill(_jets[3].pT()); // 4th jet pT
173          else if (ihist == 5 && jet_n[ithres] >= 5) _histogram(ihist, ithres)->fill(_jets[4].pT()); // 5th jet pT
174        }
175      }
176    }
177
178
179    void finalize() {
180      // Normalize to cross-section
181      const double norm = crossSection()/sumOfWeights();
182      typedef map<unsigned int, Histo1DPtr>::value_type IDtoHisto1DPtr; ///< @todo Remove when C++11 allowed
183      for (IDtoHisto1DPtr ihpair : _hMap) scale(ihpair.second, norm); ///< @todo Use normalize(ihpair.second, crossSection())
184      // Calc averages
185      for (unsigned int ihist = 0; ihist < _histLimit ; ihist++) {
186        unsigned int threshLimit = _thresholdLimit(ihist);
187        for (unsigned int ithres = 0; ithres < threshLimit; ithres++) {
188          scale(_histogram(ihist, ithres), 0.5);
189        }
190      }
191    }
192
193
194
195  private:
196
197
198    /// @name Cut helper functions
199    //@{
200
201    // Event selection functions
202    bool _ejets(unsigned int& cutBits) {
203      // 1. More than zero good electrons
204      cutBits += 1; if (_dressedelectrons.size() == 0) return false;
205      // 2. No additional electrons passing the veto selection
206      cutBits += 1 << 1; if (_vetodressedelectrons.size() > 1) return false;
207      // 3. No muons passing the veto selection
208      cutBits += 1 << 2; if (_vetodressedmuons.size() > 0) return false;
209      // 4. total neutrino pT > 30 GeV
210      cutBits += 1 << 3; if (_met_et <= 30.0*GeV) return false;
211      // 5. MTW > 35 GeV
212      cutBits += 1 << 4;
213      if (_transMass(_dressedelectrons[0].pT(), _dressedelectrons[0].phi(), _met_et, _met_phi) <= 35*GeV) return false;
214      // 6. At least one b-tagged jet
215      cutBits += 1 << 5; if (_jet_ntag < 1) return false;
216      // 7. At least three good jets
217      cutBits += 1 << 6; if (_jets.size() < 3) return false;
218      cutBits += 1 << 7;
219      return true;
220    }
221
222    bool _mujets(unsigned int& cutBits) {
223      // 1. More than zero good muons
224      cutBits += 1; if (_dressedmuons.size() == 0) return false;
225      // 2. No additional muons passing the veto selection
226      cutBits += 1 << 1; if (_vetodressedmuons.size() > 1) return false;
227      // 3. No electrons passing the veto selection
228      cutBits += 1 << 2; if (_vetodressedelectrons.size() > 0) return false;
229      // 4. total neutrino pT > 30 GeV
230      cutBits += 1 << 3; if (_met_et <= 30*GeV) return false;
231      // 5. MTW > 35 GeV
232      cutBits += 1 << 4;
233      if (_transMass(_dressedmuons[0].pT(), _dressedmuons[0].phi(), _met_et, _met_phi) <= 35*GeV) return false;
234      // 6. At least one b-tagged jet
235      cutBits += 1 << 5; if (_jet_ntag < 1) return false;
236      // 7. At least three good jets
237      cutBits += 1 << 6; if (_jets.size() < 3) return false;
238      cutBits += 1 << 7;
239      return true;
240    }
241
242    bool _additionalJetCuts() {
243      if (_jets.size() < 2) return false;
244      if (_jets[0].pT() <= 50*GeV || _jets[1].pT() <= 35*GeV) return false;
245      return true;
246    }
247
248    //@}
249
250
251    /// @name Histogram helper functions
252    //@{
253
254    unsigned int _thresholdLimit(unsigned int histId) {
255      if (histId == 0) return 4;
256      return 1;
257    }
258
259    Histo1DPtr _histogram(unsigned int histId, unsigned int thresholdId) {
260      assert(histId < _histLimit);
261      assert(thresholdId < _thresholdLimit(histId));
262
263      const unsigned int hInd = (histId == 0) ? thresholdId : (_thresholdLimit(0) + (histId-1) + thresholdId);
264      if (_hMap.find(hInd) != _hMap.end()) return _hMap[hInd];
265
266      _hMap.insert( { hInd, Histo1DPtr() } );
267      if (histId == 0) book(_hMap[hInd], thresholdId+1, 1, 1);
268      else             book(_hMap[hInd], 4+histId,      1, 1);
269      return _hMap[hInd];
270    }
271
272    //@}
273
274
275    /// @name Physics object helper functions
276    //@{
277
278    double _transMass(double ptLep, double phiLep, double met, double phiMet) {
279      return sqrt(2.0*ptLep*met*(1 - cos(phiLep-phiMet)));
280    }
281
282    unsigned int _countJets(unsigned int ithres) {
283      if (ithres > 4) assert(0);
284      double pTcut[4] = {25.,40.,60.,80.};
285      unsigned int i, jet_n = 0;
286      for (i = 0; i < _jets.size(); i++) {
287        if (_jets[i].pT() > pTcut[ithres]) jet_n++;
288      }
289      unsigned int ncutoff[4] = {8,7,6,5};
290      if (jet_n > ncutoff[ithres]) jet_n = ncutoff[ithres];
291      return jet_n;
292    }
293
294    //@}
295
296
297  private:
298
299    /// @name Objects that are used by the event selection decisions
300    //@{
301    vector<DressedLepton> _dressedelectrons;
302    vector<DressedLepton> _vetodressedelectrons;
303    vector<DressedLepton> _dressedmuons;
304    vector<DressedLepton> _vetodressedmuons;
305    Particles _neutrinos;
306    Jets _jets;
307    unsigned int _jet_ntag;
308    /// @todo Why not store the whole MET FourMomentum?
309    double _met_et, _met_phi;
310    bool _overlap;
311
312    map<unsigned int, Histo1DPtr> _hMap;
313    //unsigned int _chanLimit;
314    unsigned int _histLimit;
315    //@}
316
317  };
318
319
320
321  // Declare the class as a hook for the plugin system
322  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1304688);
323
324}