rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2017_I1609448

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