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              MSG_DEBUG("Constructing Rmiss for " << rmisslabel);
307              if (_e.find(rmisslabel) == _e.end())  book(_e[rmisslabel], rmisslabel); // book for first weight
308              // load SM prediction
309              // inject BSM component
310              YODA::Estimate1D numer = injectBSM("sr0l_"s+obs);
311              YODA::Estimate1D denom = injectBSM("sr0l_"s+auxil+"0v_"s+obs);
312              // construct Rmiss
313              if (numer != denom)  numer.rebinXTo(denom.xEdges()); // harmonise binning if need be
314              YODA::Estimate1D ratio = numer / denom; // assumes uncorrelated error breakdowns
315              // copy over old Rmiss error breakdown
316              const string path = _e[rmisslabel]->path();
317              *_e[rmisslabel] = refData(rmisslabel+"_thy_nlo"s);
318              _e[rmisslabel]->setPath(path);
319              // update central value and stats error component
320              for (auto& b : _e[rmisslabel]->bins()) {
321                const auto& injb = ratio.bin(b.index());
322                b.set(injb.val(), injb.err("stats"), "stats");
323              }
324            }
325          }
326        }
327      }
328
329
330      YODA::Estimate1D injectBSM(const string& label) {
331        YODA::Estimate1D SM = refData(label+"_thy_nlo"s);
332        for (auto& b : SM.bins())  b.scale(b.dVol());
333        return SM + _h[label]->mkEstimate("", "stats", false);
334      }
335
336      // check if jet is between tagging jets
337      bool isBetween(const Jet& probe, const Jet& boundary1, const Jet& boundary2) const {
338        double y_p = probe.rapidity();
339        double y_b1 = boundary1.rapidity();
340        double y_b2 = boundary2.rapidity();
341
342        double y_min = std::min(y_b1, y_b2);
343        double y_max = std::max(y_b1, y_b2);
344
345        return  (y_p > y_min && y_p < y_max);
346      }
347
348
349      // count number of gap jets for central jet veto
350      size_t centralJetVeto(const Jets &jets) const {
351        if (jets.size() < 2) return 0;
352        const Jet& bj1 = jets.at(0);
353        const Jet& bj2 = jets.at(1);
354
355        size_t n_between = 0;
356        // start loop at the 3rd hardest pT jet
357        for (size_t i = 2; i < jets.size(); ++i) {
358          const Jet& j = jets.at(i);
359          if (isBetween(j, bj1, bj2))  ++n_between;
360        }
361        return n_between;
362      }
363
364
365      bool hasZcandidate(const DressedLeptons& leptons) const {
366        if (leptons.size() != 2)  return false; // ask for exactly two leptons
367        if (leptons[0].pid() != -leptons[1].pid())  return false; // check for SFOS
368        double Zmass = (leptons[0].mom() + leptons[1].mom()).mass();
369        return inRange(Zmass, 66*GeV,  116*GeV); // check dilepton mass
370      }
371
372
373      double signedDeltaPhi(const Jet& j1, const Jet& j2) const {
374        double dphijj;
375        if (j1.rap() > j2.rap()) {
376          dphijj = j1.phi() - j2.phi();
377        }
378        else {
379          dphijj = j2.phi() - j1.phi();
380        }
381        return mapAngleMPiToPi(dphijj);
382      }
383
384
385      struct Categories {
386        map<string, bool> regions;
387
388        Categories () { }
389
390        void init(const string& name) { regions[name] = false; }
391
392        void trigger(const string& name) {
393          regions[name] = true;
394        }
395
396        void reset() {
397          for (auto& reg : regions) { reg.second = false; }
398        }
399      };
400
401
402    private:
403
404      // Data members like post-cuts event weight counters go here
405      /// @name Histograms
406      //@{
407      map<string, Histo1DPtr> _h;
408
409      map<string, Estimate1DPtr> _e;
410
411      Categories cats;
412
413      Cut ecal_cuts;
414
415      size_t _extra, _isBSM;
416      //@}
417
418  };
419
420  // The hook for the plugin system
421  RIVET_DECLARE_PLUGIN(ATLAS_2024_I2765017);
422}