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
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"

namespace Rivet {


  /// ATLAS pTmiss+jets cross-section ratios at 13 TeV
  class ATLAS_2017_I1609448 : public Analysis {
  public:

    /// Constructor
    ATLAS_2017_I1609448(const string name="ATLAS_2017_I1609448", 
                        const string ref_data="ATLAS_2017_I1609448") : Analysis(name) {
      setRefDataName(ref_data);
    }


    struct HistoHandler {
      Histo1DPtr histo;
      Scatter2DPtr scatter;
      unsigned int d, x, y;

      void fill(double value) {
        histo->fill(value);
      }
    };


    /// Initialize
    void init() {

      // Get options from the new option system
      _mode = 0;
      if ( getOption("LMODE") == "NU" ) _mode = 0; // using Z -> nunu channel by default
      if ( getOption("LMODE") == "MU" ) _mode = 1;
      if ( getOption("LMODE") == "EL" ) _mode = 2;

      // Prompt photons
      PromptFinalState photon_fs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 4.9);
      // Prompt electrons
      PromptFinalState el_fs(Cuts::abseta < 4.9 && Cuts::abspid == PID::ELECTRON);
      // Prompt muons
      PromptFinalState mu_fs(Cuts::abseta < 4.9 && Cuts::abspid == PID::MUON);

      // Dressed leptons
      Cut lep_cuts = Cuts::pT > 7*GeV && Cuts::abseta < 2.5;
      DressedLeptons dressed_leps(photon_fs, (_mode == 2 ? el_fs : mu_fs), 0.1, lep_cuts);
      declare(dressed_leps, "DressedLeptons");

      // In-acceptance leptons for lepton veto
      PromptFinalState veto_lep_fs(Cuts::abseta < 4.9 && (Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON));
      veto_lep_fs.acceptTauDecays();
      veto_lep_fs.acceptMuonDecays();
      DressedLeptons veto_lep(photon_fs, veto_lep_fs, 0.1, lep_cuts);
      declare(veto_lep, "VetoLeptons");

      // MET
      VetoedFinalState met_fs(Cuts::abseta > 2.5 && Cuts::abspid == PID::MUON); // veto out-of-acceptance muons
      if (_mode) met_fs.addVetoOnThisFinalState(dressed_leps);
      declare(MissingMomentum(met_fs), "MET");

      // Jet collection
      FastJets jets(FinalState(Cuts::abseta < 4.9), FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE);
      declare(jets, "Jets");

      _h["met_mono"] = bookHandler(1, 1, 2);
      _h["met_vbf" ] = bookHandler(2, 1, 2);
      _h["mjj_vbf" ] = bookHandler(3, 1, 2);
      _h["dphijj_vbf"] = bookHandler(4, 1, 2);
    }


    HistoHandler bookHandler(unsigned int id_d, unsigned int id_x, unsigned int id_y) {
      HistoHandler dummy;
      if (_mode < 2) {  // numerator mode
        const string histName = "_" + mkAxisCode(id_d, id_x, id_y);
        book(dummy.histo, histName, refData(id_d, id_x, id_y)); // hidden auxiliary output
        book(dummy.scatter, id_d, id_x, id_y - 1, true); // ratio
        dummy.d = id_d;
        dummy.x = id_x;
        dummy.y = id_y;
      } else {
        book(dummy.histo, id_d, id_x, 4); // denominator mode
      }
      return dummy;
    }


    bool isBetweenJets(const Jet& probe, const Jet& boundary1, const Jet& boundary2) {
      const double y_p = probe.rapidity();
      const double y_b1 = boundary1.rapidity();
      const double y_b2 = boundary2.rapidity();
      const double y_min = std::min(y_b1, y_b2);
      const double y_max = std::max(y_b1, y_b2);
      return (y_p > y_min && y_p < y_max);
    }


    int centralJetVeto(Jets& jets) {
      if (jets.size() < 2) return 0;
      const Jet bj1 = jets.at(0);
      const Jet bj2 = jets.at(1);

      // Start loop at the 3rd hardest pT jet
      int n_between = 0;
      for (size_t i = 2; i < jets.size(); ++i) {
        const Jet j = jets.at(i);
        if (isBetweenJets(j, bj1, bj2) && j.pT() > 25*GeV)  ++n_between;
      }
      return n_between;
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      // Require 0 (Znunu) or 2 (Zll) dressed leptons
      bool isZll = bool(_mode);
      const vector<DressedLepton> &vetoLeptons = applyProjection<DressedLeptons>(event, "VetoLeptons").dressedLeptons();
      const vector<DressedLepton> &all_leps = applyProjection<DressedLeptons>(event, "DressedLeptons").dressedLeptons();
      if (!isZll && vetoLeptons.size())    vetoEvent;
      if ( isZll && all_leps.size() != 2)  vetoEvent;

      vector<DressedLepton> leptons;
      bool pass_Zll = true;
      if (isZll) {
        // Sort dressed leptons by pT
        if (all_leps[0].pt() > all_leps[1].pt()) {
          leptons.push_back(all_leps[0]);
          leptons.push_back(all_leps[1]);
        } else {
          leptons.push_back(all_leps[1]);
          leptons.push_back(all_leps[0]);
        }
        // Leading lepton pT cut
        pass_Zll &= leptons[0].pT() > 80*GeV;
        // Opposite-charge requirement
        pass_Zll &= charge3(leptons[0]) + charge3(leptons[1]) == 0;
        // Z-mass requirement
        const double Zmass = (leptons[0].mom() + leptons[1].mom()).mass();
        pass_Zll &= (Zmass >= 66*GeV && Zmass <= 116*GeV);
      }
      if (!pass_Zll)  vetoEvent;


      // Get jets and remove those within dR = 0.5 of a dressed lepton
      Jets jets = applyProjection<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::absrap < 4.4);
      for (const DressedLepton& lep : leptons)
        ifilter_discard(jets, deltaRLess(lep, 0.5));

      const size_t njets = jets.size();
      if (!njets)  vetoEvent;
      const int njets_gap = centralJetVeto(jets);

      double jpt1 = jets[0].pT();
      double jeta1 = jets[0].eta();
      double mjj = 0., jpt2 = 0., dphijj = 0.;
      if (njets >= 2) {
        mjj = (jets[0].momentum() + jets[1].momentum()).mass();
        jpt2 = jets[1].pT();
        dphijj = deltaPhi(jets[0], jets[1]);
      }

      // MET
      Vector3 met_vec = apply<MissingMomentum>(event, "MET").vectorMPT();
      double met = met_vec.mod();

      // Cut on deltaPhi between MET and first 4 jets, but only if jet pT > 30 GeV
      bool dphi_fail = false;
      for (size_t i = 0; i < jets.size() && i < 4; ++i) {
        dphi_fail |= (deltaPhi(jets[i], met_vec) < 0.4 && jets[i].pT() > 30*GeV);
      }

      const bool pass_met_dphi = met > 200*GeV && !dphi_fail;
      const bool pass_vbf = pass_met_dphi && mjj > 200*GeV && jpt1 > 80*GeV && jpt2 > 50*GeV && njets >= 2 && !njets_gap;
      const bool pass_mono = pass_met_dphi && jpt1 > 120*GeV && fabs(jeta1) < 2.4;
      if (pass_mono)  _h["met_mono"].fill(met);
      if (pass_vbf) {
        _h["met_vbf"].fill(met/GeV);
        _h["mjj_vbf"].fill(mjj/GeV);
        _h["dphijj_vbf"].fill(dphijj);
      }
    }


    /// Normalise, scale and otherwise manipulate histograms here
    void finalize() {
      const double sf(crossSection() / femtobarn / sumOfWeights());
      for (map<string, HistoHandler>::iterator hit = _h.begin(); hit != _h.end(); ++hit) {
        scale(hit->second.histo, sf);
        if (_mode < 2)  constructRmiss(hit->second);
      }
    }


    void constructRmiss(HistoHandler& handler) {
      // Load transfer function from reference data file
      const YODA::Scatter2D& rmiss = refData(handler.d, handler.x, handler.y);
      const YODA::Scatter2D& numer = refData(handler.d, handler.x, handler.y + 1);
      const YODA::Scatter2D& denom = refData(handler.d, handler.x, handler.y + 2);
      for (size_t i = 0; i < handler.scatter->numPoints(); ++i) {
        const Point2D& r = rmiss.point(i); // SM Rmiss
        const Point2D& n = numer.point(i); // SM numerator
        const Point2D& d = denom.point(i); // SM denominator
        const HistoBin1D& b = handler.histo->bin(i); // BSM
        double bsmy;
        try {
          bsmy = b.height();
        } catch (const Exception&) { // LowStatsError or WeightError
          bsmy = 0;
        }
        double bsmey;
        try {
          bsmey = b.heightErr();
        } catch (const Exception&) { // LowStatsError or WeightError
          bsmey = 0;
        }
        // Combined numerator
        double sm_plus_bsm = n.y() + bsmy;
        // Rmiss central value
        double rmiss_y = safediv(sm_plus_bsm, d.y());
        // Ratio error (Rmiss = SM_num/SM_denom + BSM/SM_denom ~ Rmiss_SM + BSM/SM_denom
        double rmiss_p = sqrt(r.yErrPlus()*r.yErrPlus()   + safediv(bsmey*bsmey, d.y()*d.y()));
        double rmiss_m = sqrt(r.yErrMinus()*r.yErrMinus() + safediv(bsmey*bsmey, d.y()*d.y()));
        // Set new values
        Point2D& p = handler.scatter->point(i); // (SM + BSM) Rmiss
        p.setY(rmiss_y);
        p.setYErrMinus(rmiss_m);
        p.setYErrPlus(rmiss_p);
      }
    }


  protected:

    // Analysis-mode switch
    size_t _mode;

    /// Histograms
    map<string, HistoHandler> _h;

  };


  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1609448);
}

Edit this page on GitLab