rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2018_I1711223

Observation of electroweak WZjj production
Experiment: ATLAS (LHC)
Inspire ID: 1711223
Status: VALIDATED
Authors:
  • Eirini Kasimi
  • Emmanuel Sauvan
No references listed
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp -> WZ j j, diboson decays to electrons and muons, no b-quarks in the initial state

An observation of electroweak $W^\pm Z$ production in association with two jets in proton-proton collisions is presented. The data collected by the ATLAS detector at the Large Hadron Collider in 2015 and 2016 at a centre-of-mass energy of $\sqrt{s}=13$ TeV are used, corresponding to an integrated luminosity of 36.1fb$^{-1}$. Events containing three identified leptons, either electrons or muons, and two jets are selected. The electroweak production of $W^\pm Z$ bosons in association with two jets is measured with an observed significance of 5.3 standard deviations. A fiducial cross-section for electroweak production including interference effects and for a single leptonic decay mode is measured to be $\sigma_{WZjj-EW}=0.57-0.13+0.14$(stat.)$-0.06+0.07$(syst.)fb. Total and differential fiducial cross-sections of the sum of $W^\pm Zjj$ electroweak and strong productions for several kinematic observables are also measured. 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_2018_I1711223.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
 10
 11namespace Rivet {
 12
 13
 14  /// @brief Electroweak WZjj production cross section at 13 TeV
 15  class ATLAS_2018_I1711223 : public Analysis {
 16  public:
 17   
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1711223);
 20
 21    /// @name Analysis methods
 22    //@{
 23    
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27      // Get photons to dress leptons
 28      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 29
 30      // Electrons and muons in Fiducial PS
 31      PromptFinalState leptons(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
 32      leptons.acceptTauDecays(false);
 33      DressedLeptons dressedleptons(photons, leptons, 0.1, Cuts::open(), true);                                       // useDecayPhotons=true -- useJetClustering? auto-set to false?
 34      declare(dressedleptons, "DressedLeptons");
 35
 36      // Prompt neutrinos (yikes!)
 37      IdentifiedFinalState nu_id;
 38      nu_id.acceptNeutrinos();
 39      PromptFinalState neutrinos(nu_id);
 40      neutrinos.acceptTauDecays(false);
 41      declare(neutrinos, "Neutrinos");
 42      MSG_WARNING("\033[91;1mLIMITED VALIDITY - check info file for details!\033[m");
 43
 44      //Jets
 45    
 46    	// Muons
 47    	PromptFinalState bare_mu(Cuts::abspid == PID::MUON, true); // true = use muons from prompt tau decays
 48    	DressedLeptons all_dressed_mu(photons, bare_mu, 0.1, Cuts::abseta < 5.0, true);
 49
 50    	// Electrons
 51    	PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, true); // true = use electrons from prompt tau decays
 52    	DressedLeptons all_dressed_el(photons, bare_el, 0.1, Cuts::abseta < 5.0, true);
 53
 54    	//Jet forming
 55    	VetoedFinalState vfs(FinalState(Cuts::abseta < 5));
 56    	vfs.addVetoOnThisFinalState(all_dressed_el);
 57    	vfs.addVetoOnThisFinalState(all_dressed_mu);
 58          
 59    	FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
 60    	declare(jets, "Jets");
 61
 62      // Book auxiliary histograms
 63      book(_h["MTWZ"],         "_mTWZ", refData( 6, 1, 1));
 64      book(_h["sumpt"],       "_sumpT", refData( 8, 1, 1));
 65      book(_h["dphiWZ"],     "_dphiWZ", refData(10, 1, 1));
 66      book(_h["Njets_VBS"],   "_njets", refData(12, 1, 1));
 67      book(_h["mjj"],           "_mjj", refData(14, 1, 1));
 68      book(_h["dyjj"],       "_dRapjj", refData(16, 1, 1));
 69      book(_h["dphijj"],     "_dPhijj", refData(18, 1, 1));
 70      book(_h["Njets_gap"], "_gapJets", refData(20, 1, 1));
 71
 72      // book output bar charts
 73      book(_s["MTWZ"],       6, 1, 1);
 74      book(_s["sumpt"],      8, 1, 1);
 75      book(_s["dphiWZ"],    10, 1, 1);
 76      book(_s["Njets_VBS"], 12, 1, 1);
 77      book(_s["mjj"],       14, 1, 1);
 78      book(_s["dyjj"],      16, 1, 1);
 79      book(_s["dphijj"],    18, 1, 1);
 80      book(_s["Njets_gap"], 20, 1, 1);
 81      
 82    }
 83
 84
 85    void analyze(const Event& event) {
 86
 87      const Particles& dressedleptons = apply<DressedLeptons>(event, "DressedLeptons").particlesByPt();
 88      const Particles& neutrinos = apply<PromptFinalState>(event, "Neutrinos").particlesByPt();
 89      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 4.5);
 90
 91      int i, j, k;
 92      double MassZ01 = 0., MassZ02 = 0., MassZ12 = 0.;
 93      double MassW0 = 0., MassW1 = 0., MassW2 = 0.;
 94      double WeightZ1, WeightZ2, WeightZ3;
 95      double WeightW1, WeightW2, WeightW3;
 96      double M1, M2, M3;
 97      double WeightTotal1, WeightTotal2, WeightTotal3;
 98
 99      //---Fiducial PS: assign leptons to W and Z bosons using Resonant shape algorithm
100      if (dressedleptons.size() < 3 || neutrinos.size() < 1) vetoEvent;                                                              
101 
102      //--- count num of electrons and muons
103      int Nel = 0, Nmu = 0;
104      for (const Particle& l : dressedleptons) {
105        if (l.abspid() == 11)  ++Nel;
106        if (l.abspid() == 13)  ++Nmu;
107      }
108
109      int icomb=0;
110      // try Z pair of leptons 01                                                              
111    if ( (dressedleptons[0].pid() ==-(dressedleptons[1].pid()))  && (dressedleptons[2].pid()*neutrinos[0].pid()< 0) && (dressedleptons[2].abspid()==neutrinos[0].abspid()-1)) {
112        MassZ01 = (dressedleptons[0].momentum() + dressedleptons[1].momentum()).mass();
113        MassW2 = (dressedleptons[2].momentum() + neutrinos[0].momentum()).mass();
114        icomb = 1;
115      }
116
117      // try Z pair of leptons 02
118      if ( (dressedleptons[0].pid()==-(dressedleptons[2].pid()))  && (dressedleptons[1].pid()*neutrinos[0].pid()< 0) && (dressedleptons[1].abspid()==neutrinos[0].abspid()-1)) {
119        MassZ02 = (dressedleptons[0].momentum() + dressedleptons[2].momentum()).mass();
120        MassW1 = (dressedleptons[1].momentum() + neutrinos[0].momentum()).mass();
121        icomb = 2;
122      }
123      // try Z pair of leptons 12
124      if ( (dressedleptons[1].pid()==-(dressedleptons[2].pid())) && (dressedleptons[0].pid()*neutrinos[0].pid()< 0) && (dressedleptons[0].abspid()==neutrinos[0].abspid()-1)) {
125        MassZ12 = (dressedleptons[1].momentum() + dressedleptons[2].momentum()).mass();
126        MassW0 = (dressedleptons[0].momentum() + neutrinos[0].momentum()).mass();
127        icomb = 3;
128      }
129 
130      if (icomb<=0)  vetoEvent;
131
132
133      WeightZ1 = 1/(pow(MassZ01*MassZ01 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
134      WeightW1 = 1/(pow(MassW2*MassW2 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
135      WeightTotal1 = WeightZ1*WeightW1;
136      M1 = -1*WeightTotal1;
137
138      WeightZ2 = 1/(pow(MassZ02*MassZ02- MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
139      WeightW2 = 1/(pow(MassW1*MassW1- MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
140      WeightTotal2 = WeightZ2*WeightW2;
141      M2 = -1*WeightTotal2;
142
143      WeightZ3 = 1/(pow(MassZ12*MassZ12 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
144      WeightW3 = 1/(pow(MassW0*MassW0 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
145      WeightTotal3 = WeightZ3*WeightW3;
146      M3 = -1*WeightTotal3;
147
148      if( (M1 < M2 && M1 < M3) || (MassZ01 != 0 && MassW2 != 0 && MassZ02 == 0 && MassZ12 == 0) ) {
149        i = 0; j = 1; k = 2;
150      }
151      if((M2 < M1 && M2 < M3) || (MassZ02 != 0 && MassW1 != 0 && MassZ01 == 0 && MassZ12 == 0) ) {
152        i = 0; j = 2; k = 1;
153      }
154      if((M3 < M1 && M3 < M2) || (MassZ12 != 0 && MassW0 != 0 && MassZ01 == 0 && MassZ02 == 0) ) {
155        i = 1; j = 2; k = 0;
156      }
157
158      FourMomentum Zlepton1 = dressedleptons[i].momentum();
159      FourMomentum Zlepton2 = dressedleptons[j].momentum();
160      FourMomentum Wlepton  = dressedleptons[k].momentum();
161      FourMomentum Zboson   = dressedleptons[i].momentum()+dressedleptons[j].momentum();
162      FourMomentum Wboson   = dressedleptons[k].momentum()+neutrinos[0].momentum();
163
164      double cosLepNeut;
165      double Wboson_mT = 0.;
166      double norm = Wlepton.pT() * neutrinos[0].pt();
167      if(norm != 0 ) {
168        cosLepNeut = ( Wlepton.px()*neutrinos[0].px() + Wlepton.py()*neutrinos[0].py() )/norm ;
169        if (1-cosLepNeut >= 0. ) Wboson_mT = sqrt( 2 * Wlepton.pT() * neutrinos[0].pt() * (1-cosLepNeut ) );
170      }
171
172      //---- CUTS (based on Table 1 WZ: 36.1 fb-1)----//
173      if (Wlepton.pT() <= 20*GeV || Zlepton1.pT() <= 15*GeV || Zlepton2.pT() <= 15*GeV)     vetoEvent;      
174      if (Wlepton.abseta() >= 2.5 || Zlepton1.abseta() >= 2.5 || Zlepton2.abseta() >= 2.5)  vetoEvent;
175      if (fabs(Zboson.mass()/GeV - MZ_PDG) >= 10.) vetoEvent;
176      if (Wboson_mT <= 30*GeV)                     vetoEvent;
177      if (deltaR(Zlepton1, Zlepton2) <= 0.2)       vetoEvent;
178      if (deltaR(Zlepton1, Wlepton)  <= 0.3)       vetoEvent;
179      if (deltaR(Zlepton2, Wlepton)  <= 0.3)       vetoEvent;
180
181      double WZ_pt = (Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt() + neutrinos[0].pt())/GeV;
182      double WZ_px = (Zlepton1.px() + Zlepton2.px() + Wlepton.px() + neutrinos[0].px())/GeV;
183      double WZ_py = (Zlepton1.py() + Zlepton2.py() + Wlepton.py() + neutrinos[0].py())/GeV;
184      double mTWZ = sqrt( pow(WZ_pt, 2) - ( pow(WZ_px, 2) + pow(WZ_py,2) ) );
185      double sumptleptons = (Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt())/GeV;
186      double dPhiWZTruth = acos(cos(Zboson.phi()-Wboson.phi()));
187
188    
189      
190      //---- Jet CUTS----//
191      ifilter_discard(jets, [&](const Jet& j) {
192        return deltaR(j, Zlepton1) < 0.3 || deltaR(j, Zlepton2) < 0.3 || deltaR(j, Wlepton) < 0.3;
193      });
194      if (jets.size() < 2)  vetoEvent;  
195      if (jets[0].pT() < 40*GeV)  vetoEvent;  
196
197      // Selection of the second jet as the second highest pT jet and in opposite hemisphere with the fisrt jet
198      FourMomentum jet_lead = jets[0].mom();
199      FourMomentum jet_sublead;
200      bool foundVBSJetPair = false;
201      for (const Jet& jet : jets) {
202        if(jet.pT() > 40*GeV && jet.eta()*jets[0].eta() < 0.) {
203          jet_sublead = jet.mom();
204          foundVBSJetPair = true;
205          break;
206        }
207      }
208      if (!foundVBSJetPair)  vetoEvent;
209     
210      const double mJJ = (jet_lead + jet_sublead).mass()/GeV;
211      const double dphi_jj = acos(cos(jet_lead.phi() - jet_sublead.phi()));
212      const double dyjj = fabs(jet_lead.rap() - jet_sublead.rap());
213
214      //Plots in the SR
215      if (mJJ < 500*GeV) vetoEvent;
216
217      const size_t njets40 = filter_select(jets, Cuts::pT > 40*GeV).size();
218      fillWithOverflow("Njets_VBS", njets40, 5.1);
219
220      const double y_min = std::min(jet_lead.rap(), jet_sublead.rap());
221      const double y_max = std::max(jet_lead.rap(), jet_sublead.rap());
222      const size_t njetsGap = count(jets, [&](const Jet& j) { 
223        return  (j.rap() > y_min && j.rap() < y_max);
224      });
225      fillWithOverflow("Njets_gap", njetsGap, 3.1);
226
227      fillWithOverflow("MTWZ", mTWZ, 551);
228      fillWithOverflow("sumpt", sumptleptons, 501);
229      fillWithOverflow("mjj", mJJ, 2001);
230
231      _h["dphiWZ"]->fill(dPhiWZTruth);
232      _h["dyjj"]->fill(dyjj);
233      _h["dphijj"]->fill(dphi_jj);
234
235    }
236
237
238    void fillWithOverflow(const string& tag, const double value, const double overflow) {
239      if (value < overflow)  _h[tag]->fill(value);
240      else                   _h[tag]->fill(overflow);
241    }
242
243
244    void finalize() {
245
246      scale(_h, crossSectionPerEvent() / femtobarn);
247      // unfortunately, no differential cross-sections were measured in this analysis
248      for (auto &item : _h)  barchart(item.second, _s[item.first]);
249
250    }
251
252
253 //@}
254
255  private:
256
257
258    /// @name Histograms
259    //@{
260
261    map<string, Histo1DPtr> _h;
262    map<string, Scatter2DPtr> _s;
263
264    //@}
265
266    double MZ_PDG = 91.1876;
267    double MW_PDG = 80.385;
268    double GammaZ_PDG = 2.4952;
269    double GammaW_PDG = 2.085;
270
271  };
272
273  // The hook for the plugin system
274  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1711223);
275}