rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2017_I1609448

$p_\text{T}^\text{miss}$+jets cross-section ratios at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1609448
Status: VALIDATED
Authors:
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp -> MET + jets or pp -> e+ e- + jets or pp -> mu+ mu- + jets at 13 TeV

Observables sensitive to the anomalous production of events containing hadronic jets and missing momentum in the plane transverse to the proton beams at the Large Hadron Collider are presented. The observables are defined as a ratio of cross sections, for events containing jets and large missing transverse momentum to events containing jets and a pair of charged leptons from the decay of a $Z/\gamma^\ast$ boson. This definition minimises experimental and theoretical systematic uncertainties in the measurements. This ratio is measured differentially with respect to a number of kinematic properties of the hadronic system in two phase-space regions; one inclusive single-jet region and one region sensitive to vector-boson-fusion topologies. The data are found to be in agreement with the Standard Model predictions and used to constrain a variety of theoretical models for dark-matter production, including simplified models, effective field theory models, and invisible decays of the Higgs boson. The measurements use 3.2 fb${}^{-1}$ of proton-proton collision data recorded by the ATLAS experiment at a centre-of-mass energy of 13 TeV and are fully corrected for detector effects, meaning that the data can be used to constrain new-physics models beyond those shown in this paper. The reference data file comes with the measured $R^\text{miss}$ (y01), the expected $R^\text{miss}$ in the Standard Model (y02), the expected $R^\text{miss}$ numerator in the Standard Model (y03) as well as the expected $R^\text{miss}$ denominator in the Standard Model (y04). If no mode is specified, will assume routine is being run on BSM model and attempt to combine with SM prediction from ref data file. If NU is specified will assume SM $Z\to\nu\nu$ is being run on. If EL or MU is specified, will assume routine is run on SM $Z\to ee$ or $Z\to\mu\mu$ respectively. Note on reinterpretation: If a BSM signal is due to an excess of Z bosons, it should appear both numerator and denominator and so the ratio should have zero sensitivity. This will be wrongly evaluated if the SM denominator mode is used.

Source code: ATLAS_2017_I1609448.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8#include "Rivet/Projections/MissingMomentum.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// ATLAS pTmiss+jets cross-section ratios at 13 TeV
 14  class ATLAS_2017_I1609448 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1609448);
 19
 20
 21    struct HistoHandler {
 22      Histo1DPtr histo;
 23      Estimate1DPtr estimate;
 24      unsigned int d, x, y;
 25
 26      void fill(double value) {
 27        histo->fill(value);
 28      }
 29    };
 30
 31
 32    /// Initialize
 33    void init() {
 34
 35      // Get options from the new option system
 36      _mode = 0;
 37      if ( getOption("LMODE") == "NU" ) _mode = 0; // using Z -> nunu channel by default
 38      if ( getOption("LMODE") == "MU" ) _mode = 1;
 39      if ( getOption("LMODE") == "EL" ) _mode = 2;
 40
 41      // Prompt photons
 42      PromptFinalState photon_fs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 4.9);
 43      // Prompt electrons
 44      PromptFinalState el_fs(Cuts::abseta < 4.9 && Cuts::abspid == PID::ELECTRON);
 45      // Prompt muons
 46      PromptFinalState mu_fs(Cuts::abseta < 4.9 && Cuts::abspid == PID::MUON);
 47
 48      // Dressed leptons
 49      Cut lep_cuts = Cuts::pT > 7*GeV && Cuts::abseta < 2.5;
 50      LeptonFinder dressed_leps((_mode == 2 ? el_fs : mu_fs), photon_fs, 0.1, lep_cuts);
 51      declare(dressed_leps, "LeptonFinder");
 52
 53      // In-acceptance leptons for lepton veto
 54      PromptFinalState veto_lep_fs(Cuts::abseta < 4.9 && (Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON));
 55      veto_lep_fs.acceptTauDecays();
 56      veto_lep_fs.acceptMuonDecays();
 57      LeptonFinder veto_lep(veto_lep_fs, photon_fs, 0.1, lep_cuts);
 58      declare(veto_lep, "VetoLeptons");
 59
 60      // MET
 61      VetoedFinalState met_fs(Cuts::abseta > 2.5 && Cuts::abspid == PID::MUON); // veto out-of-acceptance muons
 62      if (_mode) met_fs.addVetoOnThisFinalState(dressed_leps);
 63      declare(MissingMomentum(met_fs), "MET");
 64
 65      // Jet collection
 66      FastJets jets(FinalState(Cuts::abseta < 4.9), JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 67      declare(jets, "Jets");
 68
 69      _h["met_mono"] = bookHandler(1, 1, 2);
 70      _h["met_vbf" ] = bookHandler(2, 1, 2);
 71      _h["mjj_vbf" ] = bookHandler(3, 1, 2);
 72      _h["dphijj_vbf"] = bookHandler(4, 1, 2);
 73    }
 74
 75
 76    HistoHandler bookHandler(unsigned int id_d, unsigned int id_x, unsigned int id_y) {
 77      HistoHandler dummy;
 78      if (_mode < 2) {  // numerator mode
 79        const string histName = "_" + mkAxisCode(id_d, id_x, id_y);
 80        book(dummy.histo, histName, refData(id_d, id_x, id_y)); // hidden auxiliary output
 81        book(dummy.estimate, id_d, id_x, id_y - 1); // ratio
 82        dummy.d = id_d;
 83        dummy.x = id_x;
 84        dummy.y = id_y;
 85      } else {
 86        book(dummy.histo, id_d, id_x, 4); // denominator mode
 87      }
 88      return dummy;
 89    }
 90
 91
 92    bool isBetweenJets(const Jet& probe, const Jet& boundary1, const Jet& boundary2) {
 93      const double y_p = probe.rapidity();
 94      const double y_b1 = boundary1.rapidity();
 95      const double y_b2 = boundary2.rapidity();
 96      const double y_min = std::min(y_b1, y_b2);
 97      const double y_max = std::max(y_b1, y_b2);
 98      return (y_p > y_min && y_p < y_max);
 99    }
100
101
102    int centralJetVeto(Jets& jets) {
103      if (jets.size() < 2) return 0;
104      const Jet bj1 = jets.at(0);
105      const Jet bj2 = jets.at(1);
106
107      // Start loop at the 3rd hardest pT jet
108      int n_between = 0;
109      for (size_t i = 2; i < jets.size(); ++i) {
110        const Jet j = jets.at(i);
111        if (isBetweenJets(j, bj1, bj2) && j.pT() > 25*GeV)  ++n_between;
112      }
113      return n_between;
114    }
115
116
117    /// Perform the per-event analysis
118    void analyze(const Event& event) {
119
120      // Require 0 (Znunu) or 2 (Zll) dressed leptons
121      bool isZll = bool(_mode);
122      const DressedLeptons &vetoLeptons = apply<LeptonFinder>(event, "VetoLeptons").dressedLeptons();
123      const DressedLeptons &all_leps = apply<LeptonFinder>(event, "LeptonFinder").dressedLeptons();
124      if (!isZll && vetoLeptons.size())    vetoEvent;
125      if ( isZll && all_leps.size() != 2)  vetoEvent;
126
127      DressedLeptons leptons;
128      bool pass_Zll = true;
129      if (isZll) {
130        // Sort dressed leptons by pT
131        if (all_leps[0].pt() > all_leps[1].pt()) {
132          leptons.push_back(all_leps[0]);
133          leptons.push_back(all_leps[1]);
134        } else {
135          leptons.push_back(all_leps[1]);
136          leptons.push_back(all_leps[0]);
137        }
138        // Leading lepton pT cut
139        pass_Zll &= leptons[0].pT() > 80*GeV;
140        // Opposite-charge requirement
141        pass_Zll &= charge3(leptons[0]) + charge3(leptons[1]) == 0;
142        // Z-mass requirement
143        const double Zmass = (leptons[0].mom() + leptons[1].mom()).mass();
144        pass_Zll &= (Zmass >= 66*GeV && Zmass <= 116*GeV);
145      }
146      if (!pass_Zll)  vetoEvent;
147
148
149      // Get jets and remove those within dR = 0.5 of a dressed lepton
150      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::absrap < 4.4);
151      for (const DressedLepton& lep : leptons)
152        idiscard(jets, deltaRLess(lep, 0.5));
153
154      const size_t njets = jets.size();
155      if (!njets)  vetoEvent;
156      const int njets_gap = centralJetVeto(jets);
157
158      double jpt1 = jets[0].pT();
159      double jeta1 = jets[0].eta();
160      double mjj = 0., jpt2 = 0., dphijj = 0.;
161      if (njets >= 2) {
162        mjj = (jets[0].momentum() + jets[1].momentum()).mass();
163        jpt2 = jets[1].pT();
164        dphijj = deltaPhi(jets[0], jets[1]);
165      }
166
167      // MET
168      Vector3 met_vec = apply<MissingMomentum>(event, "MET").vectorMPT();
169      double met = met_vec.mod();
170
171      // Cut on deltaPhi between MET and first 4 jets, but only if jet pT > 30 GeV
172      bool dphi_fail = false;
173      for (size_t i = 0; i < jets.size() && i < 4; ++i) {
174        dphi_fail |= (deltaPhi(jets[i], met_vec) < 0.4 && jets[i].pT() > 30*GeV);
175      }
176
177      const bool pass_met_dphi = met > 200*GeV && !dphi_fail;
178      const bool pass_vbf = pass_met_dphi && mjj > 200*GeV && jpt1 > 80*GeV && jpt2 > 50*GeV && njets >= 2 && !njets_gap;
179      const bool pass_mono = pass_met_dphi && jpt1 > 120*GeV && fabs(jeta1) < 2.4;
180      if (pass_mono)  _h["met_mono"].fill(met);
181      if (pass_vbf) {
182        _h["met_vbf"].fill(met/GeV);
183        _h["mjj_vbf"].fill(mjj/GeV);
184        _h["dphijj_vbf"].fill(dphijj);
185      }
186    }
187
188
189    /// Normalise, scale and otherwise manipulate histograms here
190    void finalize() {
191      const double sf(crossSection() / femtobarn / sumOfWeights());
192      for (auto& item : _h) {
193        scale(item.second.histo, sf);
194        if (_mode < 2)  constructRmiss(item.second);
195      }
196    }
197
198
199    void constructRmiss(HistoHandler& handler) {
200      // Load transfer function from reference data file
201      const YODA::Estimate1D& rmiss = refData(handler.d, handler.x, handler.y);
202      const YODA::Estimate1D& numer = refData(handler.d, handler.x, handler.y + 1);
203      const YODA::Estimate1D& denom = refData(handler.d, handler.x, handler.y + 2);
204      for (size_t i = 1; i < handler.estimate->numBins()+1; ++i) {
205        const auto& r = rmiss.bin(i); // SM Rmiss
206        const auto& n = numer.bin(i); // SM numerator
207        const auto& d = denom.bin(i); // SM denominator
208        const auto& b = handler.histo->bin(i); // BSM
209        double bsmy = b.sumW();
210        double bsmey = b.errW();
211        // Combined numerator
212        double sm_plus_bsm = n.val() + bsmy;
213        // Rmiss central value
214        double rmiss_y = safediv(sm_plus_bsm, d.val());
215        // Ratio error (Rmiss = SM_num/SM_denom + BSM/SM_denom ~ Rmiss_SM + BSM/SM_denom
216        double rmiss_p = sqrt(sqr(r.errPos())  + safediv(sqr(bsmey), sqr(d.val())));
217        double rmiss_m = sqrt(sqr(r.errNeg()) + safediv(sqr(bsmey), sqr(d.val())));
218        // Set new values
219        handler.estimate->bin(i).set(rmiss_y, {rmiss_m, rmiss_p});
220      }
221    }
222
223
224  protected:
225
226    // Analysis-mode switch
227    size_t _mode;
228
229    /// Histograms
230    map<string, HistoHandler> _h;
231
232  };
233
234
235  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1609448);
236}