rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2019_I1720442

Inclusive 4-lepton lineshape at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1720442
Status: VALIDATED
Authors:
  • Max Goblirsch
  • Jon Butterworth
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp+X at 13 TeV

A measurement of the four-lepton invariant mass spectrum is made with the ATLAS detector, using an integrated luminosity of 36.1 fb1 of proton-proton collisions at s=13 TeV delivered by the Large Hadron Collider. The differential cross-section is measured for events containing two same-flavour opposite-sign lepton pairs. It exhibits a rich structure, with different mass regions dominated in the Standard Model by single Z boson production, Higgs boson production, and Z boson pair production, and non-negligible interference effects at high invariant masses. The measurement is compared with state-of-the-art Standard Model calculations, which are found to be consistent with the data. These calculations are used to interpret the data in terms of ggZZ4 and Z4 subprocesses, and to place constraints on a possible contribution from physics beyond the Standard Model.

Source code: ATLAS_2019_I1720442.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// ATLAS 4-lepton lineshape at 13 TeV
 11  class ATLAS_2019_I1720442 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1720442);
 16
 17    void init() {
 18
 19      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 20      PromptFinalState elecs(Cuts::abspid == PID::ELECTRON);
 21      PromptFinalState muons(Cuts::abspid == PID::MUON);
 22
 23      // Selection
 24      Cut el_fid_sel = (Cuts::abseta < 2.47) && (Cuts::pT > 7*GeV);
 25      Cut mu_fid_sel = (Cuts::abseta < 2.7) && (Cuts::pT > 5*GeV);
 26
 27      LeptonFinder dressed_elecs(elecs, photons, 0.005, el_fid_sel);
 28      declare(dressed_elecs, "elecs");
 29
 30      LeptonFinder dressed_muons(muons, photons, 0.005, mu_fid_sel);
 31      declare(dressed_muons, "muons");
 32
 33      // Book histos
 34      book(_h["m4l_inclusive"], 1,1,1);
 35
 36      book(_h["m4l_ptslice1"], 2,1,1);
 37      book(_h["m4l_ptslice2"], 3,1,1);
 38      book(_h["m4l_ptslice3"], 4,1,1);
 39      book(_h["m4l_ptslice4"], 5,1,1);
 40
 41      book(_h["m4l_rapidityslice1"], 6,1,1);
 42      book(_h["m4l_rapidityslice2"], 7,1,1);
 43      book(_h["m4l_rapidityslice3"], 8,1,1);
 44      book(_h["m4l_rapidityslice4"], 9,1,1);
 45
 46      book(_h["m4l_4mu"], 12,1,1);
 47      book(_h["m4l_4e"], 13,1,1);
 48      book(_h["m4l_2e2mu"], 14,1,1);
 49    }
 50
 51
 52    /// @brief Generic dilepton candidate
 53    /// @todo Move into the Rivet core?
 54    struct Dilepton : public ParticlePair {
 55      Dilepton() { }
 56      Dilepton(ParticlePair _particlepair) : ParticlePair(_particlepair) {
 57        assert(first.abspid() == second.abspid());
 58      }
 59      FourMomentum mom() const { return first.momentum() + second.momentum(); }
 60      operator FourMomentum() const { return mom(); }
 61      static bool cmppT(const Dilepton& lx, const Dilepton& rx) { return lx.mom().pT() < rx.mom().pT(); }
 62      int flavour() const { return first.abspid(); }
 63      double pTl1() const { return first.pT(); }
 64      double pTl2() const { return second.pT(); }
 65    };
 66
 67
 68    struct Quadruplet {
 69      Quadruplet (Dilepton z1, Dilepton z2): _z1(z1), _z2(z2) { }
 70      enum class FlavCombi { mm=0, ee, me, em, undefined };
 71      FourMomentum mom() const { return _z1.mom() + _z2.mom(); }
 72      Dilepton getZ1() const { return _z1; }
 73      Dilepton getZ2() const { return _z2; }
 74      Dilepton _z1, _z2;
 75      FlavCombi type() const {
 76        if (     _z1.flavour() == 13 && _z2.flavour() == 13) { return FlavCombi::mm; }
 77        else if (_z1.flavour() == 11 && _z2.flavour() == 11) { return FlavCombi::ee; }
 78        else if (_z1.flavour() == 13 && _z2.flavour() == 11) { return FlavCombi::me; }
 79        else if (_z1.flavour() == 11 && _z2.flavour() == 13) { return FlavCombi::em; }
 80        else  return FlavCombi::undefined;
 81      }
 82    };
 83
 84
 85    vector<Quadruplet> getBestQuads(Particles& particles) {
 86      // H->ZZ->4l pairing
 87      // - Two same flavor opposite charged leptons
 88      // - Ambiguities in pairing are resolved by choosing the combination
 89      //     that results in the smaller value of |mll - mZ| for each pair successively
 90      vector<Quadruplet> quads {};
 91
 92      size_t n_parts = particles.size();
 93      if (n_parts < 4)  return quads;
 94
 95      // STEP 1: find SFOS pairs
 96      vector<Dilepton> SFOS;
 97      for (size_t i = 0; i < n_parts; ++i) {
 98        for (size_t j = 0; j < i; ++j) {
 99          if (particles[i].pid() == -particles[j].pid()) {
100            // sort such that the negative lepton is listed first
101            if (particles[i].pid() > 0)  SFOS.push_back(Dilepton(make_pair(particles[i], particles[j])));
102            else                         SFOS.push_back(Dilepton(make_pair(particles[j], particles[i])));
103          }
104        }
105      }
106      if (SFOS.size() < 2)  return quads;
107
108      // now we sort the SFOS pairs
109      std::sort(SFOS.begin(), SFOS.end(), [](const Dilepton& p1, const Dilepton& p2) {
110          return fabs(p1.mom().mass() - Z_mass) < fabs(p2.mom().mass() - Z_mass);
111        });
112
113      // Form all possible quadruplets passing the pT cuts
114      for (size_t k = 0; k < SFOS.size(); ++k) {
115        for (size_t l = k+1; l < SFOS.size(); ++l) {
116          if (deltaR(SFOS[k].first.mom(),  SFOS[l].first.mom())  < 1e-13)  continue;
117          if (deltaR(SFOS[k].first.mom(),  SFOS[l].second.mom()) < 1e-13)  continue;
118          if (deltaR(SFOS[k].second.mom(), SFOS[l].first.mom())  < 1e-13)  continue;
119          if (deltaR(SFOS[k].second.mom(), SFOS[l].second.mom()) < 1e-13)  continue;
120
121          vector<double> lep_pt { SFOS[k].pTl1(), SFOS[k].pTl2(), SFOS[l].pTl1(), SFOS[l].pTl2() };
122          std::sort(lep_pt.begin(), lep_pt.end(), std::greater<double>());
123          if (!(lep_pt[0] > 20*GeV && lep_pt[1] > 15*GeV && lep_pt[2] > 10*GeV)) continue;
124          quads.push_back( Quadruplet(SFOS[k], SFOS[l]) );
125        }
126      }
127      return quads;
128    }
129
130
131    bool passMassCuts(const Quadruplet& theQuad){
132      const vector<double> cuts_m34{ 5*GeV, 5*GeV, 12*GeV, 12*GeV, 50*GeV };
133      const vector<double> cuts_m4l{ 0, 100*GeV, 110*GeV, 140*GeV, 190*GeV };
134
135      double m4l = theQuad.mom().mass();
136      double mZ1 = theQuad.getZ1().mom().mass();
137      double mZ2 = theQuad.getZ2().mom().mass();
138
139      // Invariant-mass requirements
140      double cutval = cuts_m34.back();
141      for (size_t k = 0; k < cuts_m34.size(); ++k) {
142        if (cuts_m4l[k] > m4l) {
143          cutval = cuts_m34[k-1] + (cuts_m34[k] - cuts_m34[k-1])/(cuts_m4l[k] - cuts_m4l[k-1]) * (m4l - cuts_m4l[k-1]);
144          break;
145        }
146      }
147      return inRange(mZ1, 50*GeV, 106*GeV) && inRange(mZ2, cutval, 115*GeV);
148    }
149
150
151    bool pass_dRll(const Quadruplet& theQuad) {
152      const double dR_min_same = 0.1;
153      const double dR_min_opp = 0.2;
154      double dr_min_cross = dR_min_opp;
155      if (theQuad.getZ1().flavour() == theQuad.getZ2().flavour()) {
156        dr_min_cross = dR_min_same;
157      }
158      return !((deltaR(theQuad.getZ1().first,  theQuad.getZ1().second) < dR_min_same)  ||
159               (deltaR(theQuad.getZ2().first,  theQuad.getZ2().second) < dR_min_same)  ||
160               (deltaR(theQuad.getZ1().first,  theQuad.getZ2().first)  < dr_min_cross) ||
161               (deltaR(theQuad.getZ1().first,  theQuad.getZ2().second) < dr_min_cross) ||
162               (deltaR(theQuad.getZ1().second, theQuad.getZ2().first)  < dr_min_cross) ||
163               (deltaR(theQuad.getZ1().second, theQuad.getZ2().second) < dr_min_cross));
164    }
165
166
167    bool pass_Jpsi(const Quadruplet& theQuad){
168      Particles all_leps { theQuad.getZ1().first, theQuad.getZ1().second, theQuad.getZ2().first, theQuad.getZ2().second };
169      for (const Particle& lep1 : all_leps) {
170        for (const Particle& lep2 : all_leps) {
171          if (lep1.pid() == -lep2.pid() && (lep1.mom() + lep2.mom()).mass() < 5*GeV) return false;
172        }
173      }
174      return true;
175    }
176
177
178    // Handle 3 further CF stages - m12/34, dRmin, jpsi veto
179    bool passSelection (const Quadruplet& theQuad){
180      return passMassCuts(theQuad) && pass_dRll(theQuad) && pass_Jpsi(theQuad);
181    }
182
183
184    // Do the analysis
185    void analyze(const Event& event) {
186
187      //preselection of leptons for ZZ-> llll final state
188      Particles dressed_leptons;
189      for (auto lep : apply<LeptonFinder>(event, "muons").dressedLeptons()) { dressed_leptons.push_back(lep); }
190      for (auto lep : apply<LeptonFinder>(event, "elecs").dressedLeptons()) { dressed_leptons.push_back(lep); }
191
192      auto foundDressed = getBestQuads(dressed_leptons);
193      // if we don't find any quad, we can stop here
194      if (foundDressed.empty())  vetoEvent;
195
196      bool pass = passSelection(foundDressed[0]);
197      if (pass) {
198        double m4l = foundDressed[0].mom().mass()/GeV;
199        double pt4l = foundDressed[0].mom().pT()/GeV;
200        double y4l = foundDressed[0].mom().absrap();
201        Quadruplet::FlavCombi flavour = foundDressed[0].type();
202        _h["m4l_inclusive"]->fill(m4l);
203        if (     pt4l <  20.)  _h["m4l_ptslice1"]->fill(m4l);
204        else if (pt4l <  50.)  _h["m4l_ptslice2"]->fill(m4l);
205        else if (pt4l < 100.)  _h["m4l_ptslice3"]->fill(m4l);
206        else if (pt4l < 600.)  _h["m4l_ptslice4"]->fill(m4l);
207
208        if (     y4l < 0.4)  _h["m4l_rapidityslice1"]->fill(m4l);
209        else if (y4l < 0.8)  _h["m4l_rapidityslice2"]->fill(m4l);
210        else if (y4l < 1.2)  _h["m4l_rapidityslice3"]->fill(m4l);
211        else if (y4l < 2.5)  _h["m4l_rapidityslice4"]->fill(m4l);
212
213        if (     flavour == Quadruplet::FlavCombi::mm) _h["m4l_4mu"]->fill(m4l);
214        else if (flavour == Quadruplet::FlavCombi::ee) _h["m4l_4e"]->fill(m4l);
215        else if (flavour == Quadruplet::FlavCombi::me || flavour == Quadruplet::FlavCombi::em) {
216          _h["m4l_2e2mu"]->fill(m4l);
217        }
218      }
219
220    }
221
222
223    /// Finalize
224    void finalize() {
225      const double sf = crossSection() / femtobarn / sumOfWeights();
226      scale(_h, sf);
227    }
228
229
230  private:
231
232    map<string, Histo1DPtr> _h;
233    static constexpr double Z_mass = 91.1876;
234
235  };
236
237
238  RIVET_DECLARE_PLUGIN(ATLAS_2019_I1720442);
239
240}