rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2024_I2765017

pTmiss+jets cross-sections and ratios at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 2765017
Status: VALIDATED
Authors:
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp -> 0/1/2 l + jet(s) (l = e, mu, nu) at 13 TeV

Measurements of inclusive, differential cross-sections for the production of events with missing transverse momentum in association with jets in proton-proton collisions at $\sqrt{s} =$ 13 TeV are presented. The measurements are made with the ATLAS detector using an integrated luminosity of 140 fb$^{-1}$ and include measurements of dijet distributions in a region in which vector-boson fusion processes are enhanced. They are unfolded to correct for detector resolution and efficiency within the fiducial acceptance, and are designed to allow robust comparisons with a wide range of theoretical predictions. A measurement of differential cross sections for the $Z\to\nu\nu$ process is made. The measurements are generally well-described by Standard Model predictions except for the dijet invariant mass distribution. Auxiliary measurements of the hadronic system recoiling against isolated leptons, and photons, are also made in the same phase space. Ratios between the measured distributions are then derived, to take advantage of cancellations in modelling effects and some of the major systematic uncertainties. These measurements are sensitive to new phenomena, and provide a mechanism to easily set constraints on phenomenological models. To illustrate the robustness of the approach, these ratios are compared with two common Dark Matter models, where the constraints derived from the measurement are comparable to those set by dedicated detector-level searches.

Source code: ATLAS_2024_I2765017.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#include "Rivet/Projections/InvisibleFinalState.hh"
 10#include "Rivet/Projections/VisibleFinalState.hh"
 11
 12namespace Rivet {
 13
 14  /// @brief pTmiss+jets cross-sections and ratios at 13 TeV
 15  class ATLAS_2024_I2765017 : public Analysis {
 16
 17    public:
 18
 19      /// Constructor
 20      RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2024_I2765017);
 21
 22    public:
 23
 24      // initialise
 25      void init() {
 26
 27        _extra = getOption("MODE", "FIDUCIAL") == "EXTRAPOLATED";
 28        _isBSM = getOption("TYPE", "BSM") != "SM";
 29
 30        ecal_cuts = Cuts::abseta < 2.47 && !(Cuts::abseta > 1.37 && Cuts::abseta < 1.52);
 31
 32        // prompt photons
 33        PromptFinalState photons(Cuts::abspid == PID::PHOTON && Cuts::abseta < 4.9);
 34        declare(photons, "Photons");
 35
 36        // prompt (bare) leptons
 37        Cut lep_inc_cuts = Cuts::abseta < 4.9 && (Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
 38        PromptFinalState leps(lep_inc_cuts, TauDecaysAs::PROMPT);
 39
 40        // dressed leptons
 41        Cut el_fid_cuts  = (ecal_cuts && Cuts::abspid == PID::ELECTRON);
 42        Cut mu_fid_cuts  = (Cuts::abseta < 2.5 && Cuts::abspid == PID::MUON);
 43        Cut lep_fid_cuts = (Cuts::pT > 7*GeV) && (mu_fid_cuts || el_fid_cuts);
 44        LeptonFinder dressed_leps(leps, photons, 0.1, lep_fid_cuts);
 45        declare(dressed_leps, "DressedLeptons");
 46
 47        // jet collection
 48        VetoedFinalState jet_fs(FinalState(Cuts::abseta < 4.9));
 49        if (_extra)  jet_fs.addVetoOnThisFinalState(dressed_leps);
 50        FastJets jets(jet_fs, JetAlg::ANTIKT, 0.4,
 51                      _extra? JetMuons::DECAY : JetMuons::NONE,
 52                      _extra? JetInvisibles::DECAY : JetInvisibles::NONE);
 53        declare(jets, "Jets");
 54
 55        // MET
 56        if (_extra) declare(InvisibleFinalState(OnlyPrompt::YES), "promptNeutrinos");
 57        else {
 58          FinalState met_fs(!(Cuts::abspid == PID::MUON && (Cuts::abseta > 2.5 || Cuts::pT < 7*GeV)));
 59          declare(MissingMomentum(met_fs), "actualMET");
 60        }
 61
 62        // Calorimeter particles for photon isolation
 63        VetoedFinalState calo_fs(VisibleFinalState(Cuts::abspid != PID::MUON));
 64        declare(calo_fs, "CaloParticles");
 65
 66        // Jets for UE subtraction with jet-area method
 67        if (_extra) {
 68          FastJets ktjets(FinalState(), JetAlg::KT, 0.5, JetMuons::NONE, JetInvisibles::NONE);
 69          ktjets.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec(0.9)));
 70          declare(ktjets, "kTjets");
 71        }
 72
 73        // Book Histograms
 74        cats.init("sr0l");
 75        cats.init("sr0l_cr1ie_0v"); cats.init("sr0l_cr1im_0v");
 76        cats.init("sr0l_cr2ie_0v"); cats.init("sr0l_cr2im_0v");
 77        if (_extra) cats.init("sr0l_cr1ig_0v");
 78
 79        for (const auto& cat : cats.regions) {
 80          const string mod = cat.first + "_";
 81
 82          // HepData-based booking
 83          const string phsp1 = "met_mono";
 84          const string phsp2 = "met_vbf";
 85          const string phsp3 = "mjj_vbf";
 86          const string phsp4 = "dphijj_vbf";
 87          const string phsp4dbp = "dphij_vbf"; // Silly typo :(
 88          const string suff = _extra? "_dbp" : "";
 89          // book the ones in the paper from yoda file
 90          if (mod=="sr0l_") {
 91            book(_h[mod + phsp1], mod + phsp1 + suff);
 92            book(_h[mod + phsp2], mod + phsp2 + suff);
 93            book(_h[mod + phsp3], mod + phsp3 + suff);
 94            book(_h[mod + phsp4], mod + phsp4 + suff);
 95          } else if (mod=="sr0l_cr2im_0v_") {
 96            book(_h[mod + phsp1], mod + phsp1 + suff);
 97            book(_h[mod + phsp2], mod + phsp2 + suff);
 98            book(_h[mod + phsp3], mod + phsp3 + suff);
 99            book(_h[mod + phsp4], mod + phsp4 + suff);
100          } else if (mod=="sr0l_cr1im_0v_") {
101            book(_h[mod + phsp1], mod + phsp1 + suff);
102            book(_h[mod + phsp2], mod + phsp2 + suff);
103            book(_h[mod + phsp3], mod + phsp3 + suff);
104            book(_h[mod + phsp4], mod + phsp4 + suff);
105          } else if (mod=="sr0l_cr2ie_0v_") {
106            book(_h[mod + phsp1], mod + phsp1 + suff);
107            book(_h[mod + phsp2], mod + phsp2 + suff);
108            book(_h[mod + phsp3], mod + phsp3 + suff);
109            book(_h[mod + phsp4], mod + (_extra? phsp4dbp : phsp4) + suff);
110          } else if (mod=="sr0l_cr1ie_0v_") {
111            book(_h[mod + phsp1], mod + phsp1 + suff);
112            book(_h[mod + phsp2], mod + phsp2 + suff);
113            book(_h[mod + phsp3], mod + phsp3 + suff);
114            book(_h[mod + phsp4], mod + phsp4 + suff);
115          }
116          else if (_extra && mod=="sr0l_cr1ig_0v_") {
117            book(_h[mod + phsp1], mod + phsp1 + suff);
118            book(_h[mod + phsp2], mod + phsp2 + suff);
119            book(_h[mod + phsp3], mod + phsp3 + suff);
120            book(_h[mod + phsp4], mod + phsp4 + suff);
121          }
122
123        }
124      } // end of initialize
125
126      /// Perform the per-event analysis
127      void analyze(const Event& event) {
128
129        cats.reset();
130
131        // get prompt photons
132        const Particles &photons = apply<PromptFinalState>(event, "Photons").particlesByPt(Cuts::abseta < 2.47 && Cuts::pT > 7*GeV);
133
134        // get dressed leptons
135        DressedLeptons leptons = apply<LeptonFinder>(event, "DressedLeptons").dressedLeptons(cmpMomByPt);
136        // additional lepton veto
137        if (leptons.size() > 2)  vetoEvent;
138
139        // veto on *prompt* hadronic tau candidates
140        for (const auto& jet : apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && ecal_cuts)) {
141          for (const auto& p : jet.particles()) {
142            if (p.fromHadronicTau(true))  vetoEvent; // true = only consider prompt taus
143          }
144        }
145
146        // get anti-kT 0.4 jets
147        Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4);
148
149        // remove jets within dR = 0.2 cone of a dressed lepton
150        idiscardIfAnyDeltaRLess(jets, leptons, 0.2);
151        // remove dressed leptons within dR = 0.4 cone of a jet
152        idiscardIfAnyDeltaRLess(leptons, jets, 0.4);
153        // remove jets within dR = 0.2 of a prompt photon
154        if (_extra && leptons.empty() && photons.size()) {
155          idiscardIfAnyDeltaRLess(jets, photons, 0.2);
156        }
157
158        const size_t njets = jets.size();
159        const size_t njets_gap = centralJetVeto(jets);
160
161        // calculate MET
162        Vector3 METvec;
163        if (_extra) {
164          METvec = sum(apply<InvisibleFinalState>(event, "promptNeutrinos").particles(), Kin::pTvec, METvec).setZ(0);
165        }
166        else {
167          METvec = apply<MissingMomentum>(event, "actualMET").vectorMissingPt();
168        }
169        const double actualMET = METvec.mod()/GeV; // actual pTmiss
170        Vector3 pMETvec = sum(leptons, Kin::pTvec, METvec).setZ(0); // pseudo-MET with 'invisible' leptons
171        if (_extra && leptons.empty() && photons.size()) {
172          pMETvec = (photons[0].pTvec() + pMETvec).setZ(0); // pseudo-MET with 'invisible' photon
173        }
174        const double ptmiss = pMETvec.mod()/GeV;
175
176        // lepton-MET kinematics
177        const double lep_pt = leptons.size()? leptons[0].pT()/GeV : 0.;
178        const bool hasZ = hasZcandidate(leptons);
179        double dphi_metl = -1., mT_l = 0.;
180        if (leptons.size()) {
181          dphi_metl = deltaPhi(leptons[0], METvec); // actual MET here since SR1 and SR2 have visible leps
182          mT_l = sqrt( 2 * lep_pt * actualMET * (1 - cos(dphi_metl) ) );
183          if (leptons.size() > 1) {
184            FourMomentum Z = leptons[0].mom() + leptons[1].mom();
185          }
186        }
187
188        // photon kinematics
189        const double photon_pt = photons.size()? photons[0].pT()/GeV : 0.;
190
191        // photon isolation
192        bool is_iso_gamma = true;
193        if (_extra) {
194          const vector<double> eta_bins = {0.0, 1.5, 3.0};
195          vector<double> rho(eta_bins.size()-1, 0.0);
196          FastJets kTjets = apply<FastJets>(event, "kTjets");
197          for (size_t ieta = 0; ieta < eta_bins.size()-1; ++ieta) {
198            fastjet::Selector fjselector(fastjet::SelectorAbsRapRange(eta_bins[ieta], eta_bins[ieta+1]));
199            double sigma, area;
200            kTjets.clusterSeqArea()->get_median_rho_and_sigma(fjselector, true, rho[ieta], sigma, area);
201          }
202          const double isoRCone = 0.4;
203          for (const Particle& photon : photons) {
204            // Compute calo isolation via particles within a cone around the photon
205            FourMomentum mom_in_EtCone;
206            for (const Particle& p : apply<VetoedFinalState>(event, "CaloParticles").particles()) {
207              if (deltaR(photon, p) > isoRCone)  continue; // reject if not in cone
208              mom_in_EtCone += p; // sum momentum
209            }
210            mom_in_EtCone -= photon; // subtract core photon
211            // UE subtraction energy
212            const double UEpT = M_PI * sqr(isoRCone) * rho[binIndex(photon.abseta(), eta_bins)];
213            // Use photon if energy in isolation cone is low enough
214            if (mom_in_EtCone.Et()/GeV - UEpT > 0.044*photon.mom().pT()/GeV + 2.45) is_iso_gamma = false;
215            break;
216          }
217        }
218
219        // jet kinematics
220        double jpt1 = 0, jeta1 = 0., jpt2 = 0.;
221        double mjj = 0., drapjj = 0., dphijj = 0.;
222        //double dphi_metj = 0., mT_j = 0.;
223        if (njets) {
224          jpt1 = jets[0].pT()/GeV;
225          jeta1 = jets[0].eta();
226          if (njets >= 2) {
227            mjj = (jets[0].mom() + jets[1].mom()).mass()/GeV;
228            jpt2 = jets[1].pT()/GeV;
229            drapjj = deltaRap(jets[0], jets[1]);
230            dphijj = signedDeltaPhi(jets[0], jets[1]);
231          }
232        }
233
234        // jet-MET balance (check dPhi between MET and first 4 jets)
235        bool fail_dphi_jet_MET = false;
236        for (size_t i = 0; i < jets.size() && i < 4; ++i) {
237          fail_dphi_jet_MET |= (deltaPhi(jets[i], pMETvec) < 0.4);
238        }
239
240        // start categorising
241        if (leptons.empty() && ptmiss > 200. && !fail_dphi_jet_MET) { // 0-lep SRs
242          cats.trigger("sr0l");
243        }
244
245        if (leptons.size() == 1 && ptmiss > 200. && !fail_dphi_jet_MET) { // CR for 0-lep SR using W
246          // higher leading electron pT cut due to trigger, cut on actual MET to suppress fakes
247          if (leptons[0].abspid() == PID::ELECTRON && lep_pt > 30. && actualMET > 60. && inRange(mT_l, 30., 100.)) {
248            cats.trigger("sr0l_cr1ie_0v");
249          }
250          else if (leptons[0].abspid() == PID::MUON) {
251            cats.trigger("sr0l_cr1im_0v");
252          }
253        }
254        if (hasZ && ptmiss > 200. && lep_pt > 80. && !fail_dphi_jet_MET) { // CR for 0-lep SR using Z
255          if (leptons[0].abspid() == PID::ELECTRON) {
256            cats.trigger("sr0l_cr2ie_0v");
257          }
258          else if (leptons[0].abspid() == PID::MUON) {
259            cats.trigger("sr0l_cr2im_0v");
260          }
261        }
262        if (_extra && leptons.empty() && ptmiss > 200. && photons.size() && photon_pt > 160.) {
263          if (!fail_dphi_jet_MET && is_iso_gamma) { // CR for 0-lep SR using gamma
264            cats.trigger("sr0l_cr1ig_0v");
265          }
266        }
267
268        // identify jet topology
269        const bool pass_monojet = jpt1 > 120*GeV && fabs(jeta1) < 2.4;
270        const bool pass_vbf = njets >= 2 && mjj > 200*GeV && jpt1 > 80. && jpt2 > 50. && fabs(drapjj) > 1.; // CJV applied below
271
272        // fill histograms for all categories
273        for (const auto& cat : cats.regions) {
274
275          // check if category was triggered
276          if (!cat.second)  continue;
277          // construct prefix for histogram name
278          const string mod = cat.first + "_";
279
280          if (pass_monojet) {
281            _h[mod + "met_mono"]->fill(ptmiss);
282          } // end of pass monojet
283
284          if (pass_vbf) {
285            if (!njets_gap) { // gap-jet veto
286              // This is the VBF region!
287              _h[mod + "met_vbf"]->fill(ptmiss);
288              _h[mod + "mjj_vbf"]->fill(mjj);
289              _h[mod + "dphijj_vbf"]->fill(dphijj/pi);
290            }
291          } // end of pass VBF
292
293        } // end of loop over categories
294
295      }// end of analyze
296
297
298      /// Normalise, scale and otherwise manipulate histograms here
299      void finalize() {
300        scale(_h, crossSection() / sumOfWeights() / femtobarn);
301
302        if (_isBSM && !_extra) {
303          for (const string& auxil : vector<string>{"cr1ie_"s, "cr1im_"s, "cr2ie_"s, "cr2im_"s}) {
304            for (const string& obs : vector<string>{"met_mono"s, "met_vbf"s, "mjj_vbf"s, "dphijj_vbf"s}) {
305              const string rmisslabel("rmiss_"s+auxil+obs);
306              const string statslabel("stats");
307              MSG_DEBUG("Constructing Rmiss for " << rmisslabel);
308              if (_e.find(rmisslabel) == _e.end())  book(_e[rmisslabel], rmisslabel); // book for first weight
309              // load SM prediction
310              const string numlabel("sr0l_"s+obs);
311              const string denlabel("sr0l_"s+auxil+"0v_"s+obs);
312              const YODA::Estimate1D& numSM = refData(numlabel+"_thy_nlo"s);
313              const YODA::Estimate1D& denSM = refData(denlabel+"_thy_nlo"s);
314              const YODA::Estimate1D& rmiss = refData(rmisslabel+"_thy_nlo"s);
315              // inject BSM component
316              YODA::Estimate1D numer = numSM + _h[numlabel]->mkEstimate("", statslabel);
317              YODA::Estimate1D denom = denSM + _h[denlabel]->mkEstimate("", statslabel);
318              // construct Rmiss
319              if (numer != denom)  numer.rebinXTo(denom.xEdges()); // harmonise binning if need be
320              YODA::Estimate1D ratio = numer / denom; // assumes uncorrelated error breakdowns
321              // copy over old Rmiss error breakdown
322              const string path = _e[rmisslabel]->path();
323              *_e[rmisslabel] = rmiss; _e[rmisslabel]->setPath(path);
324              // update central value and stats error component
325              for (auto& b : _e[rmisslabel]->bins()) {
326                const auto& injb = ratio.bin(b.index());
327                b.set(injb.val(), injb.err(statslabel), statslabel);
328              }
329            }
330          }
331        }
332      }
333
334
335      // check if jet is between tagging jets
336      bool isBetween(const Jet& probe, const Jet& boundary1, const Jet& boundary2) const {
337        double y_p = probe.rapidity();
338        double y_b1 = boundary1.rapidity();
339        double y_b2 = boundary2.rapidity();
340
341        double y_min = std::min(y_b1, y_b2);
342        double y_max = std::max(y_b1, y_b2);
343
344        return  (y_p > y_min && y_p < y_max);
345      }
346
347
348      // count number of gap jets for central jet veto
349      size_t centralJetVeto(const Jets &jets) const {
350        if (jets.size() < 2) return 0;
351        const Jet& bj1 = jets.at(0);
352        const Jet& bj2 = jets.at(1);
353
354        size_t n_between = 0;
355        // start loop at the 3rd hardest pT jet
356        for (size_t i = 2; i < jets.size(); ++i) {
357          const Jet& j = jets.at(i);
358          if (isBetween(j, bj1, bj2))  ++n_between;
359        }
360        return n_between;
361      }
362
363
364      bool hasZcandidate(const DressedLeptons& leptons) const {
365        if (leptons.size() != 2)  return false; // ask for exactly two leptons
366        if (leptons[0].pid() != -leptons[1].pid())  return false; // check for SFOS
367        double Zmass = (leptons[0].mom() + leptons[1].mom()).mass();
368        return inRange(Zmass, 66*GeV,  116*GeV); // check dilepton mass
369      }
370
371
372      double signedDeltaPhi(const Jet& j1, const Jet& j2) const {
373        double dphijj;
374        if (j1.rap() > j2.rap()) {
375          dphijj = j1.phi() - j2.phi();
376        }
377        else {
378          dphijj = j2.phi() - j1.phi();
379        }
380        return mapAngleMPiToPi(dphijj);
381      }
382
383
384      struct Categories {
385        map<string, bool> regions;
386
387        Categories () { }
388
389        void init(const string& name) { regions[name] = false; }
390
391        void trigger(const string& name) {
392          regions[name] = true;
393        }
394
395        void reset() {
396          for (auto& reg : regions) { reg.second = false; }
397        }
398      };
399
400
401    private:
402
403      // Data members like post-cuts event weight counters go here
404      /// @name Histograms
405      //@{
406      map<string, Histo1DPtr> _h;
407
408      map<string, Estimate1DPtr> _e;
409
410      Categories cats;
411
412      Cut ecal_cuts;
413
414      size_t _extra, _isBSM;
415      //@}
416
417  };
418
419  // The hook for the plugin system
420  RIVET_DECLARE_PLUGIN(ATLAS_2024_I2765017);
421}