rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2021_I1849535

Inclusive 4-lepton cross sections at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1849535
Status: VALIDATED
Authors:
  • Max Goblirsch
  • Jon Butterworth
  • Louie Dartmoor Corpe
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • $p p \to \ell \ell \ell \ell + X$

Measurements of four-lepton differential and integrated fiducial cross-sections in events with two same-flavour, opposite-charge electron or muon pairs are presented. The data correspond to 139 fb$^{-1}$ of $\sqrt{s}=13$ TeV proton-proton collisions, collected by the ATLAS detector during Run 2 of the Large Hadron Collider (2015--2018). The final state has contributions from a number of interesting Standard Model processes that dominate in different four-lepton invariant mass regions, including single $Z$ boson production, Higgs boson production and on-shell $ZZ$ production, with a complex mix of interference terms, and possible contributions from physics beyond the Standard Model. The differential cross-sections include the four-lepton invariant mass inclusively, in slices of other kinematic variables, and in different lepton flavour categories. Also measured are dilepton invariant masses, transverse momenta, and angular correlation variables, in four regions of four-lepton invariant mass, each dominated by different processes. The measurements are corrected for detector effects and are compared with state-of-the-art Standard Model calculations, which are found to be consistent with the data. The $Z\to 4\ell$ branching fraction is extracted, giving a value of $\left(4.41 \pm 0.30\right) \times 10^{-6}$. Constraints on effective field theory parameters and a model based on a spontaneously broken $B-L$ gauge symmetry are also evaluated. Further reinterpretations can be performed with the provided information.

Source code: ATLAS_2021_I1849535.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7#include "Rivet/Projections/DressedLeptons.hh"
  8
  9namespace Rivet {
 10
 11  /// @name M4lLineshape analysis
 12  class ATLAS_2021_I1849535 : public Analysis {
 13    public:
 14
 15      /// Constructor
 16      RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2021_I1849535);
 17
 18      void init() {
 19
 20        // Selection
 21        Cut el_fid_sel = (Cuts::abseta < 2.47) && (Cuts::pT > 7*GeV); 
 22        Cut mu_fid_sel = (Cuts::abseta < 2.7) && (Cuts::pT > 5*GeV);
 23
 24        PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 25        PromptFinalState elecs(Cuts::abspid == PID::ELECTRON);
 26        PromptFinalState muons(Cuts::abspid == PID::MUON && mu_fid_sel);
 27        elecs.acceptTauDecays(true);
 28        muons.acceptTauDecays(true);
 29
 30
 31        // Final state including all charged particles
 32        declare(ChargedFinalState(), "CFS");
 33
 34        DressedLeptons dressed_elecs(photons, elecs, 0.1, el_fid_sel, false);
 35        declare(dressed_elecs, "elecs");
 36
 37        declare(muons, "muons");
 38
 39        // Book histos
 40        book(_h["m4l_paper"],      1,1,1);
 41        book(_h["m4l_4mu_paper"],  2,1,1);
 42        book(_h["m4l_4e_paper"],   3,1,1);
 43        book(_h["m4l_2e2mu_paper"],4,1,1);
 44
 45        book(_h["mZ1_Z_paper"],5,1,1);  
 46        book(_h["mZ1_H_paper"],6,1,1);
 47        book(_h["mZ1_offshell_paper"],7,1,1);
 48        book(_h["mZ1_ZZ_paper"],8,1,1);
 49
 50        book(_h["mZ2_Z_paper"],9,1,1);
 51        book(_h["mZ2_H_paper"],10,1,1);
 52        book(_h["mZ2_offshell_paper"],11,1,1);
 53        book(_h["mZ2_ZZ_paper"],12,1,1);
 54
 55        book(_h["ptZ1_Z_paper"],13,1,1);
 56        book(_h["ptZ1_H_paper"],14,1,1);
 57        book(_h["ptZ1_offshell_paper"],15,1,1); 
 58        book(_h["ptZ1_ZZ_paper"],16,1,1);
 59
 60        book(_h["ptZ2_Z_paper"],17,1,1);
 61        book(_h["ptZ2_H_paper"],18,1,1);
 62        book(_h["ptZ2_offshell_paper"],19,1,1);
 63        book(_h["ptZ2_ZZ_paper"],20,1,1);
 64
 65        book(_h["costhetastar1_Z_paper"],21,1,1);
 66        book(_h["costhetastar1_H_paper"],22,1,1);
 67        book(_h["costhetastar1_offshell_paper"],23,1,1);
 68        book(_h["costhetastar1_ZZ_paper"],24,1,1);
 69
 70        book(_h["costhetastar2_Z_paper"],25,1,1);
 71        book(_h["costhetastar2_H_paper"],26,1,1);
 72        book(_h["costhetastar2_offshell_paper"],27,1,1);
 73        book(_h["costhetastar2_ZZ_paper"],28,1,1);
 74
 75        book(_h["dy_Z1Z2_Z_paper"]  ,29,1,1);
 76        book(_h["dy_Z1Z2_H_paper"]  ,30,1,1);
 77        book(_h["dy_Z1Z2_offshell_paper"],31,1,1);
 78        book(_h["dy_Z1Z2_ZZ_paper"]  ,32,1,1);
 79
 80        book(_h["dphi_Z1Z2_Z_paper"]  ,33,1,1);
 81        book(_h["dphi_Z1Z2_H_paper"]  ,34,1,1);
 82        book(_h["dphi_Z1Z2_offshell_paper"],35,1,1);
 83        book(_h["dphi_Z1Z2_ZZ_paper"]  ,36,1,1);
 84
 85        book(_h["dphi_l1l2_Z_paper"]   ,37,1,1);
 86        book(_h["dphi_l1l2_H_paper"]  ,38,1,1);
 87        book(_h["dphi_l1l2_offshell_paper"],39,1,1);
 88        book(_h["dphi_l1l2_ZZ_paper"]  ,40,1,1);
 89
 90        book(_h["m4l_ptslice1_paper"],41,1,1);
 91        book(_h["m4l_ptslice2_paper"],42,1,1);
 92        book(_h["m4l_ptslice3_paper"],43,1,1);
 93        book(_h["m4l_ptslice4_paper"],44,1,1);
 94        book(_h["m4l_ptslice5_paper"],45,1,1);
 95
 96        book(_h["m4l_yslice1_paper"],46,1,1);
 97        book(_h["m4l_yslice2_paper"],47,1,1);
 98        book(_h["m4l_yslice3_paper"],48,1,1);
 99        book(_h["m4l_yslice4_paper"],49,1,1);
100        book(_h["m4l_yslice5_paper"],50,1,1);
101
102      }
103
104      /// Generic dilepton candidate
105      struct Dilepton : public ParticlePair {
106        Dilepton() { }
107        Dilepton(ParticlePair _particlepair) : ParticlePair(_particlepair) {
108          assert(first.abspid() == second.abspid());
109        }
110        FourMomentum mom() const { return first.momentum() + second.momentum(); }
111        operator FourMomentum() const { return mom(); }
112        static bool cmppT(const Dilepton& lx, const Dilepton& rx) { return lx.mom().pT() < rx.mom().pT(); }
113        int flavour() const { return first.abspid(); }
114        double pTl1() const { return first.pT(); }
115        double pTl2() const { return second.pT(); }
116      };
117
118      struct Quadruplet {
119        Quadruplet (Dilepton z1, Dilepton z2): _z1(z1), _z2(z2) { }
120        enum class FlavCombi { mm=0, ee, me, em, undefined };
121        FourMomentum mom() const { return _z1.mom() + _z2.mom(); }
122        Dilepton getZ1() const { return _z1; }
123        Dilepton getZ2() const { return _z2; }
124        Dilepton _z1, _z2;
125        FlavCombi type() const {
126          if (     _z1.flavour() == 13 && _z2.flavour() == 13) { return FlavCombi::mm; }
127          else if (_z1.flavour() == 11 && _z2.flavour() == 11) { return FlavCombi::ee; } 
128          else if (_z1.flavour() == 13 && _z2.flavour() == 11) { return FlavCombi::me; }
129          else if (_z1.flavour() == 11 && _z2.flavour() == 13) { return FlavCombi::em; }
130          else  return FlavCombi::undefined;
131        }
132      };
133      bool passesTruthIsolation(Quadruplet quad, const Particles charged_tracks, Particles& truthLeptons ){
134        bool pass =true;
135        Particles leps;
136        leps.push_back(quad._z1.first);
137        leps.push_back(quad._z2.first);
138        leps.push_back(quad._z1.second);
139        leps.push_back(quad._z2.second);
140        for (auto &lep : leps){
141          double pTinCone = -lep.pT();
142          for (const Particle& track : charged_tracks) {
143            if (deltaR(lep.momentum(), track.momentum()) < 0.3)
144              pTinCone += track.pT();
145          }
146          for (const Particle& tlep: truthLeptons) {
147            float dR= deltaR(lep.momentum(),  tlep.momentum()); 
148            if ( dR>0 && dR < 0.3)
149              pTinCone -= tlep.pT();
150          }
151          if (pTinCone > 0.16* lep.pT()){
152            pass=false;
153          }
154        }
155        return pass;
156      }
157
158      std::vector<Quadruplet> getBestQuads(Particles& particles, bool drcut = true) {
159        // H->ZZ->4l pairing
160        // - Two same flavor opposite charged leptons
161        // - Ambiguities in pairing are resolved by choosing the combination
162        //     that results in the smaller value of |mll - mZ| for each pair successively
163        std::vector<Quadruplet> quads {};
164
165        size_t n_parts = particles.size();
166        if (n_parts < 4)  return quads; 
167
168        // STEP 1: find SFOS pairs 
169        std::vector<Dilepton> SFOS;
170        for (size_t i = 0; i < n_parts; ++i) {
171          for (size_t j = 0; j < i; ++j) {
172            if (particles[i].pid() == -particles[j].pid()) {
173              // sort such that the negative lepton is listed first
174              Dilepton sfos;
175              if (particles[i].pid() > 0)  sfos = make_pair(particles[i], particles[j]);
176              else                         sfos = make_pair(particles[j], particles[i]);
177
178              if (sfos.mom().mass() > _ll_mass*GeV && (!drcut || deltaR(particles[i],particles[j]) > _dRll)) SFOS.push_back(sfos);
179            }
180          }
181        }
182        if (SFOS.size() < 2)  return quads;
183
184        // now we sort the SFOS pairs
185        std::sort(SFOS.begin(), SFOS.end(), [](const Dilepton& p1, const Dilepton& p2) {
186            return fabs(p1.mom().mass() - Z_mass) < fabs(p2.mom().mass() - Z_mass);
187            });
188
189        //form all possible quadruplets, passing the pt cuts, the dR cuts and the mll cuts
190        for (size_t k = 0; k < SFOS.size(); ++k) {
191          for (size_t l = k+1; l < SFOS.size(); ++l) {
192            if(drcut) {
193              if (deltaR(SFOS[k].first.mom(),  SFOS[l].first.mom())  < _dRll)  continue;
194              if (deltaR(SFOS[k].first.mom(),  SFOS[l].second.mom()) < _dRll)  continue;
195              if (deltaR(SFOS[k].second.mom(), SFOS[l].first.mom())  < _dRll)  continue;
196              if (deltaR(SFOS[k].second.mom(), SFOS[l].second.mom()) < _dRll)  continue;
197            }
198            if ( (SFOS[k].first.pid()   == -SFOS[l].first.pid())  && ((SFOS[k].first.mom()  + SFOS[l].first.mom()).mass()  < _ll_mass*GeV)) continue;
199            if ( (SFOS[k].first.pid()   == -SFOS[l].second.pid()) && ((SFOS[k].first.mom()  + SFOS[l].second.mom()).mass() < _ll_mass*GeV)) continue;
200            if ( (SFOS[k].second.pid()  == -SFOS[l].first.pid())  && ((SFOS[k].second.mom() + SFOS[l].first.mom()).mass()  < _ll_mass*GeV)) continue;
201            if ( (SFOS[k].second.pid()  == -SFOS[l].second.pid()) && ((SFOS[k].second.mom() + SFOS[l].second.mom()).mass() < _ll_mass*GeV)) continue;
202
203            //think technically this should happen before quad formation now and with all leptons not just those in quad so commenting out
204            //std::vector<double> lep_pt { SFOS[k].pTl1(), SFOS[k].pTl2(), SFOS[l].pTl1(), SFOS[l].pTl2() };
205            //std::sort(lep_pt.begin(), lep_pt.end(), std::greater<double>()); 
206            //if (!(lep_pt[0] > _pt_lep1*GeV && lep_pt[1] > _pt_lep2*GeV && lep_pt[2] > _pt_lep3*GeV)) continue;
207            quads.push_back( Quadruplet(SFOS[k], SFOS[l]) );
208          }
209        }
210        return quads;
211      }
212
213      bool passPtLeptons(const Particles& particles) {
214        size_t n_parts = particles.size();
215        if (n_parts < 4)  return false;
216        // cut on pT of leptons
217        return ( particles[0].mom().pt() > _pt_lep1*GeV && particles[1].mom().pt() > _pt_lep2*GeV && particles[2].mom().pt() > _pt_lep3*GeV) ;
218      }
219
220
221      // Do the analysis
222      void analyze(const Event& event) {
223
224        const Particles charged_tracks    = apply<ChargedFinalState>(event, "CFS").particles();
225
226
227        //preselection of leptons for ZZ-> llll final state
228        Particles dressed_leptons;
229        for (auto lep : apply<FinalState>(event, "muons").particles()) { dressed_leptons.push_back(lep); }
230        for (auto lep : apply<DressedLeptons>(event, "elecs").dressedLeptons()) { dressed_leptons.push_back(lep); }
231
232
233
234        // sort to put highest pT first
235        std::sort(dressed_leptons.begin(), dressed_leptons.end(), [](const Particle& l1, const Particle& l2) {
236            return l1.pt() > l2.pt();
237            });
238
239        auto foundDressedNoDrll = getBestQuads(dressed_leptons,false);
240
241        //now doing pt cut before quad formation so also apply this here
242        if (!passPtLeptons(dressed_leptons)) vetoEvent;
243
244        auto foundDressed = getBestQuads(dressed_leptons);
245        // if we don't find any quad, we can stop here 
246        if (foundDressed.empty())  vetoEvent;
247        if (!passesTruthIsolation(foundDressed[0], charged_tracks, dressed_leptons)) vetoEvent;
248
249        double m4l = foundDressed[0].mom().mass()/GeV;
250        double pt4l = foundDressed[0].mom().pT()/GeV;
251        double y4l = foundDressed[0].mom().absrap();
252        double mZ1 = foundDressed[0].getZ1().mom().mass()/GeV;
253        double mZ2 = foundDressed[0].getZ2().mom().mass()/GeV;
254        double ptZ1 = foundDressed[0].getZ1().mom().pT()/GeV;
255        double ptZ2 = foundDressed[0].getZ2().mom().pT()/GeV;
256        double dy_Z1Z2 = fabs(foundDressed[0].getZ1().mom().rapidity() - foundDressed[0].getZ2().mom().rapidity());
257        double dphi_Z1Z2 = deltaPhi(foundDressed[0].getZ1().mom(),foundDressed[0].getZ2().mom());
258        double dphi_l1l2 = deltaPhi(dressed_leptons[0].mom(),dressed_leptons[1].mom());
259
260        _h["m4l_paper"]->fill(m4l);
261        if (     pt4l <  10.)	  _h["m4l_ptslice1_paper"]->fill(m4l);
262        else if (pt4l <  20.)	  _h["m4l_ptslice2_paper"]->fill(m4l);
263        else if (pt4l < 50.)	  _h["m4l_ptslice3_paper"]->fill(m4l);
264        else if (pt4l < 100.)	  _h["m4l_ptslice4_paper"]->fill(m4l);
265        else if (pt4l < 600.)	  _h["m4l_ptslice5_paper"]->fill(m4l); 
266
267        if (y4l < 0.3)	  _h["m4l_yslice1_paper"]->fill(m4l);
268        else if (y4l < 0.6)	  _h["m4l_yslice2_paper"]->fill(m4l);
269        else if (y4l < 0.9)	  _h["m4l_yslice3_paper"]->fill(m4l);
270        else if (y4l < 1.2)	  _h["m4l_yslice4_paper"]->fill(m4l);
271        else if (y4l < 2.5)	  _h["m4l_yslice5_paper"]->fill(m4l);
272
273
274        Quadruplet::FlavCombi flavour = foundDressed[0].type();
275        if (     flavour == Quadruplet::FlavCombi::mm) { _h["m4l_4mu_paper"]->fill(m4l); }
276        else if (flavour == Quadruplet::FlavCombi::ee) { _h["m4l_4e_paper"]->fill(m4l); }
277        else if (flavour == Quadruplet::FlavCombi::me || flavour == Quadruplet::FlavCombi::em) {
278          _h["m4l_2e2mu_paper"]->fill(m4l);
279        }
280
281        // polarization variables
282        // Get four-momentum of the first lepton pair
283        const FourMomentum pcom = foundDressed.at(0).getZ1().mom();
284        const Vector3 betacom = pcom.betaVec();
285        const Vector3 unitboostvec = betacom.unit();
286        const LorentzTransform comboost = LorentzTransform::mkFrameTransformFromBeta(betacom);
287        // Get four-momentum of the negative lepton w.r.t. the first lepton pair
288        const FourMomentum p1com = comboost.transform(foundDressed.at(0).getZ1().first.mom());
289        float costhetastar1 = cos(p1com.p3().angle(unitboostvec));
290
291        // Get four-momentum of the second lepton pair
292        const FourMomentum pcom2 = foundDressed.at(0).getZ2().mom();
293        const Vector3 betacom2 = pcom2.betaVec();
294        const Vector3 unitboostvec2 = betacom2.unit();
295        const LorentzTransform comboost2 = LorentzTransform::mkFrameTransformFromBeta(betacom2);
296        // Get four-momentum of the negative lepton w.r.t. the second lepton pair
297        const FourMomentum p2com = comboost2.transform(foundDressed.at(0).getZ2().first.mom());
298        float  costhetastar2 = cos(p2com.p3().angle(unitboostvec2));
299
300        //fill m4l binned variables
301        if (60 < m4l && m4l < 100.) {
302          _h["mZ1_Z_paper"]->fill(mZ1);
303          _h["mZ2_Z_paper"]->fill(mZ2);
304          _h["ptZ1_Z_paper"]->fill(ptZ1);
305          _h["ptZ2_Z_paper"]->fill(ptZ2);
306          _h["dy_Z1Z2_Z_paper"]->fill(dy_Z1Z2);
307          _h["dphi_Z1Z2_Z_paper"]->fill(dphi_Z1Z2);
308          _h["dphi_l1l2_Z_paper"]->fill(dphi_l1l2);
309          _h["costhetastar1_Z_paper"]->fill(costhetastar1 );
310          _h["costhetastar2_Z_paper"]->fill(costhetastar2 );
311        }
312        else if(120 < m4l && m4l < 130 ){
313          _h["mZ1_H_paper"]->fill(mZ1);
314          _h["mZ2_H_paper"]->fill(mZ2);
315          _h["ptZ1_H_paper"]->fill(ptZ1);
316          _h["ptZ2_H_paper"]->fill(ptZ2);
317          _h["dy_Z1Z2_H_paper"]->fill(dy_Z1Z2);
318          _h["dphi_Z1Z2_H_paper"]->fill(dphi_Z1Z2);
319          _h["dphi_l1l2_H_paper"]->fill(dphi_l1l2);
320          _h["costhetastar1_H_paper"]->fill(costhetastar1 );
321          _h["costhetastar2_H_paper"]->fill(costhetastar2 );
322        }
323        else if(180 < m4l && m4l < 2000){
324          _h["mZ1_ZZ_paper"]->fill(mZ1);
325          _h["mZ2_ZZ_paper"]->fill(mZ2);
326          _h["ptZ1_ZZ_paper"]->fill(ptZ1);
327          _h["ptZ2_ZZ_paper"]->fill(ptZ2);
328          _h["dy_Z1Z2_ZZ_paper"]->fill(dy_Z1Z2);
329          _h["dphi_Z1Z2_ZZ_paper"]->fill(dphi_Z1Z2);
330          _h["dphi_l1l2_ZZ_paper"]->fill(dphi_l1l2);
331          _h["costhetastar1_ZZ_paper"]->fill(costhetastar1 );
332          _h["costhetastar2_ZZ_paper"]->fill(costhetastar2 );
333        }
334        else{
335          _h["mZ1_offshell_paper"]->fill(mZ1);
336          _h["mZ2_offshell_paper"]->fill(mZ2);
337          _h["ptZ1_offshell_paper"]->fill(ptZ1);
338          _h["ptZ2_offshell_paper"]->fill(ptZ2);
339          _h["dy_Z1Z2_offshell_paper"]->fill(dy_Z1Z2);
340          _h["dphi_Z1Z2_offshell_paper"]->fill(dphi_Z1Z2);
341          _h["dphi_l1l2_offshell_paper"]->fill(dphi_l1l2);
342          _h["costhetastar1_offshell_paper"]->fill(costhetastar1 );
343          _h["costhetastar2_offshell_paper"]->fill(costhetastar2 );
344        }
345
346      }//end analysis
347
348      /// Finalize
349      void finalize() {
350        const double sf = crossSection() / femtobarn / sumOfWeights();
351        for (auto hist : _h) { scale(hist.second, sf); }
352      }
353
354    private:
355
356      map<string, Histo1DPtr> _h;
357      static constexpr double Z_mass = 91.1876;
358      static constexpr float _pt_lep1 = 20.;
359      static constexpr float _pt_lep2 = 10.;
360      static constexpr float _pt_lep3 = 0.;
361      static constexpr float _ll_mass = 5.;
362      static constexpr float _dRll = 0.05;
363
364  };  // end class ATLAS_2021_I1849535
365
366  RIVET_DECLARE_PLUGIN(ATLAS_2021_I1849535);
367}  // end namespace rivet