rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2019_I1720438

Measurement of the $WZ$ production cross section at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1720438
Status: UNVALIDATED
Authors:
  • Eirini Kasimi
  • Emmanuel Sauvan
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp -> WZ + X, diboson decays to electrons or muons, data only for single channel

The production of $W^\pm Z$ events in proton--proton collisions at a centre-of-mass energy of 13 TeV is measured with the ATLAS detector at the LHC. The collected data correspond to an integrated luminosity of 36.1 fb${}^{-1}$ . The $W^\pm Z$ candidates are reconstructed using leptonic decays of the gauge bosons into electrons or muons. The measured inclusive cross section in the detector fiducial region for leptonic decay modes is $\sigma(W^\pm Z\to\ell^\prime\nu\ell\ell \text{fid.}) = 63.7\pm 3.2$(stat.)$\pm 1.0$(sys.)$\pm 1.4$(lumi.) fb. In comparison, the next-to-leading-order Standard Model prediction is 61.5-1.3+1.4 fb. Cross sections for $W^+Z$ and $W^-Z$ production and their ratio are presented as well as differentialcross sections for several kinematic observable. Users should note that explicit matching of lepton flavour between individual SM neutrinos and charged leptons is used in this analysis routine, to match the MC-based correction to the fiducial region applied in the paper. The data are therefore only valid under the assumption of the Standard Model and cannot be used for BSM reinterpretation. Uses SM neutrino-lepton flavour matching and a resonant shape algorithm assuming the Standard Model, to match the MC-based correction to the fiducial region applied in the paper. This routine is therefore only valid under the assumption of the Standard Model and cannot be used for BSM reinterpretation

Source code: ATLAS_2019_I1720438.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/DressedLeptons.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief WZ production cross-section at 13 TeV
 14  class ATLAS_2019_I1720438 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1720438);
 19
 20    /// @name Analysis methods
 21    //@{
 22
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 28
 29      // Electrons and muons in Fiducial PS
 30      PromptFinalState leptons(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
 31      leptons.acceptTauDecays(false);
 32      DressedLeptons dressedleptons(photons, leptons, 0.1, Cuts::open(), true);
 33      declare(dressedleptons, "DressedLeptons");
 34
 35      // Prompt neutrinos (yikes!)
 36      IdentifiedFinalState nu_id;
 37      nu_id.acceptNeutrinos();
 38      PromptFinalState neutrinos(nu_id);
 39      neutrinos.acceptTauDecays(false);
 40      declare(neutrinos, "Neutrinos");
 41      MSG_WARNING("\033[91;1mLIMITED VALIDITY - check info file for details!\033[m");
 42
 43      // Muons
 44      PromptFinalState bare_mu(Cuts::abspid == PID::MUON, true); // true = use muons from prompt tau decays
 45      DressedLeptons all_dressed_mu(photons, bare_mu, 0.1, Cuts::abseta < 2.5, true);
 46
 47      // Electrons
 48      PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, true); // true = use electrons from prompt tau decays
 49      DressedLeptons all_dressed_el(photons, bare_el, 0.1, Cuts::abseta < 2.5, true);
 50
 51      //Jet forming
 52      VetoedFinalState vfs(FinalState(Cuts::abseta < 5.0));
 53      vfs.addVetoOnThisFinalState(all_dressed_el);
 54      vfs.addVetoOnThisFinalState(all_dressed_mu);
 55	      
 56      FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
 57      declare(jets, "Jets");
 58
 59      // Book auxiliary histograms
 60      book(_h["pTZ"],    "_pTZ",    refData( 8, 1, 1));
 61      book(_h["pTW"],    "_pTW",    refData(10, 1, 1));
 62      book(_h["mTWZ"],   "_mTWZ",   refData(12, 1, 1));
 63      book(_h["dPhiWZ"], "_dPhiWZ", refData(14, 1, 1));
 64      book(_h["pTv"],    "_pTV",    refData(16, 1, 1));
 65      book(_h["dRapWZ"], "_drapWZ", refData(18, 1, 1));
 66      book(_h["Njets"],  "_njets",  refData(20, 1, 1)); 
 67      book(_h["Mjj"],    "_mjj",    refData(22, 1, 1)); 
 68
 69      // book output bar charts
 70      book(_s["pTZ"],     8, 1, 1);
 71      book(_s["pTW"],    10, 1, 1);
 72      book(_s["mTWZ"],   12, 1, 1);
 73      book(_s["dPhiWZ"], 14, 1, 1);
 74      book(_s["pTv"],    16, 1, 1);
 75      book(_s["dRapWZ"], 18, 1, 1);
 76      book(_s["Njets"],  20, 1, 1); 
 77      book(_s["Mjj"],    22, 1, 1); 
 78
 79    }
 80
 81
 82    void analyze(const Event& event) {
 83
 84      const Particles& dressedleptons = apply<DressedLeptons>(event, "DressedLeptons").particlesByPt();
 85      const Particles& neutrinos = apply<PromptFinalState>(event, "Neutrinos").particlesByPt();
 86      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 4.5);                                                                           
 87      int i, j, k;
 88      double MassZ01 = 0., MassZ02 = 0., MassZ12 = 0.;
 89      double MassW0 = 0., MassW1 = 0., MassW2 = 0.;
 90      double WeightZ1, WeightZ2, WeightZ3;
 91      double WeightW1, WeightW2, WeightW3;
 92      double M1, M2, M3;
 93      double WeightTotal1, WeightTotal2, WeightTotal3;
 94
 95      //---Fiducial PS: assign leptons to W and Z bosons using Resonant shape algorithm
 96       if (dressedleptons.size() < 3 || neutrinos.size() < 1) vetoEvent;                                                              
 97
 98      //--- count num of electrons and muons
 99      int Nel = 0, Nmu = 0;
100      for (const Particle& l : dressedleptons) {
101        if (l.abspid() == 11)  ++Nel;
102        if (l.abspid() == 13)  ++Nmu;
103      }
104
105      int icomb=0;
106      // try Z pair of leptons 01                                                              
107    if ( (dressedleptons[0].pid() ==-(dressedleptons[1].pid()))  && (dressedleptons[2].pid()*neutrinos[0].pid()< 0) && (dressedleptons[2].abspid()==neutrinos[0].abspid()-1)) {
108        MassZ01 = (dressedleptons[0].momentum() + dressedleptons[1].momentum()).mass();
109        MassW2 = (dressedleptons[2].momentum() + neutrinos[0].momentum()).mass();
110        icomb = 1;
111      }
112
113      // try Z pair of leptons 02
114      if ( (dressedleptons[0].pid()==-(dressedleptons[2].pid()))  && (dressedleptons[1].pid()*neutrinos[0].pid()< 0) && (dressedleptons[1].abspid()==neutrinos[0].abspid()-1)) {
115        MassZ02 = (dressedleptons[0].momentum() + dressedleptons[2].momentum()).mass();
116        MassW1 = (dressedleptons[1].momentum() + neutrinos[0].momentum()).mass();
117        icomb = 2;
118      }
119      // try Z pair of leptons 12
120      if ( (dressedleptons[1].pid()==-(dressedleptons[2].pid())) && (dressedleptons[0].pid()*neutrinos[0].pid()< 0) && (dressedleptons[0].abspid()==neutrinos[0].abspid()-1)) {
121        MassZ12 = (dressedleptons[1].momentum() + dressedleptons[2].momentum()).mass();
122        MassW0 = (dressedleptons[0].momentum() + neutrinos[0].momentum()).mass();
123        icomb = 3;
124      }
125 
126      if (icomb<=0)  vetoEvent;
127
128
129      WeightZ1 = 1/(pow(MassZ01*MassZ01 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
130      WeightW1 = 1/(pow(MassW2*MassW2 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
131      WeightTotal1 = WeightZ1*WeightW1;
132      M1 = -1*WeightTotal1;
133
134      WeightZ2 = 1/(pow(MassZ02*MassZ02- MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
135      WeightW2 = 1/(pow(MassW1*MassW1- MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
136      WeightTotal2 = WeightZ2*WeightW2;
137      M2 = -1*WeightTotal2;
138
139      WeightZ3 = 1/(pow(MassZ12*MassZ12 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
140      WeightW3 = 1/(pow(MassW0*MassW0 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
141      WeightTotal3 = WeightZ3*WeightW3;
142      M3 = -1*WeightTotal3;
143
144      if( (M1 < M2 && M1 < M3) || (MassZ01 != 0 && MassW2 != 0 && MassZ02 == 0 && MassZ12 == 0) ) {
145        i = 0; j = 1; k = 2;
146      }
147      if((M2 < M1 && M2 < M3) || (MassZ02 != 0 && MassW1 != 0 && MassZ01 == 0 && MassZ12 == 0) ) {
148        i = 0; j = 2; k = 1;
149      }
150      if((M3 < M1 && M3 < M2) || (MassZ12 != 0 && MassW0 != 0 && MassZ01 == 0 && MassZ02 == 0) ) {
151        i = 1; j = 2; k = 0;
152      }
153
154      FourMomentum Zlepton1 = dressedleptons[i].mom();
155      FourMomentum Zlepton2 = dressedleptons[j].mom();
156      FourMomentum Wlepton  = dressedleptons[k].mom();
157      FourMomentum Zboson   = dressedleptons[i].mom()+dressedleptons[j].mom();
158      FourMomentum Wboson   = dressedleptons[k].mom()+neutrinos[0].mom();
159
160     double cosLepNeut;
161     double Wboson_mT = 0;
162   	 double norm = Wlepton.pT() * neutrinos[0].pt();
163	   if (norm != 0) {
164	     cosLepNeut = ( Wlepton.px()*neutrinos[0].px() + Wlepton.py()*neutrinos[0].py() )/norm;
165	     if ( 1-cosLepNeut >= 0 )  Wboson_mT = sqrt( 2 * Wlepton.pT() * neutrinos[0].pt() * (1-cosLepNeut ) );
166	   }
167
168      //---- CUTS (based on Table 1 WZ: 36.1 fb-1)----//
169      if (Wlepton.pT() <= 20*GeV || Zlepton1.pT() <= 15*GeV || Zlepton2.pT() <= 15*GeV)  vetoEvent;      
170      if (Wlepton.abseta() >= 2.5 || Zlepton1.abseta() >= 2.5 || Zlepton2.abseta() >= 2.5)  vetoEvent;
171      if (fabs(Zboson.mass()/GeV - MZ_PDG) >= 10.) vetoEvent;
172      if (Wboson_mT <= 30*GeV)                     vetoEvent;
173      if (deltaR(Zlepton1, Zlepton2) <= 0.2)       vetoEvent;
174      if (deltaR(Zlepton1, Wlepton)  <= 0.3)       vetoEvent;
175      if (deltaR(Zlepton2, Wlepton)  <= 0.3)       vetoEvent;
176
177      double pTZ = Zboson.pT()/GeV;
178      double WZ_pt = (Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt() + neutrinos[0].pt())/GeV;
179      double WZ_px = (Zlepton1.px() + Zlepton2.px() + Wlepton.px() + neutrinos[0].px())/GeV;
180      double WZ_py = (Zlepton1.py() + Zlepton2.py() + Wlepton.py() + neutrinos[0].py())/GeV;
181      double mTWZ = sqrt( pow(WZ_pt, 2) - ( pow(WZ_px, 2) + pow(WZ_py,2) ));
182
183      double dPhiWZTruth = acos(cos(Zboson.phi()-Wboson.phi()));
184      double pTW = Wboson.pT()/GeV;
185      double pTv = neutrinos[0].pT()/GeV;
186      double AbsDeltay = fabs(Zboson.rapidity()-Wlepton.rapidity());
187
188      ifilter_discard(jets, [&](const Jet& j) {
189        return deltaR(j, Zlepton1) < 0.3 || deltaR(j, Zlepton2) < 0.3 || deltaR(j, Wlepton) < 0.3;
190      });
191
192      size_t njets = jets.size()>5? 5 : jets.size();
193      _h["Njets"]->fill(njets);
194
195      if (njets > 1) {
196        double mjj = (jets[0].mom() + jets[1].mom()).mass()/GeV;
197        if (mjj > 800.)  mjj = 800.;
198        _h["Mjj"]->fill(mjj);
199      }
200
201      if (pTZ > 220.) pTZ = 220.;
202      _h["pTZ"]->fill(pTZ);
203
204      if (pTW > 220.)  pTW = 220.;
205      _h["pTW"]->fill(pTW);
206
207      if (mTWZ > 600.)  mTWZ = 600.;
208      _h["mTWZ"]->fill(mTWZ);
209
210      _h["dPhiWZ"]->fill(dPhiWZTruth);
211
212      if (pTv > 90.)  pTv = 90.;
213      _h["pTv"]->fill(pTv);
214
215      _h["dRapWZ"]->fill(AbsDeltay);
216
217    }
218
219
220    void finalize() {
221
222      scale(_h, 0.25 * crossSectionPerEvent() / femtobarn); // data values are for _single_ lepton channel
223      // unfortunately, no differential cross-sections were measured in this analysis
224      for (auto &item : _h)  barchart(item.second, _s[item.first]);
225
226    }
227
228    //@}
229
230
231  private:
232
233
234    /// @name Histograms
235    //@{
236    
237    map<string, Histo1DPtr> _h;
238    map<string, Scatter2DPtr> _s;
239
240    //@}
241
242    double MZ_PDG = 91.1876;
243    double MW_PDG = 80.385;
244    double GammaZ_PDG = 2.4952;
245    double GammaW_PDG = 2.085;
246
247  };
248
249  // The hook for the plugin system
250  RIVET_DECLARE_PLUGIN(ATLAS_2019_I1720438);
251}