Rivet analyses referenceATLAS_2018_I1711223Observation of electroweak WZjj productionExperiment: ATLAS (LHC) Inspire ID: 1711223 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
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/LeptonFinder.hh"
8#include "Rivet/Projections/VetoedFinalState.hh"
9
10namespace Rivet {
11
12
13 /// @brief Electroweak WZjj production cross section at 13 TeV
14 class ATLAS_2018_I1711223 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1711223);
19
20 /// @name Analysis methods
21 /// @{
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25
26 // Get photons to dress leptons
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 LeptonFinder dressedleptons(leptons, photons, 0.1);
33 declare(dressedleptons, "LeptonFinder");
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 //Jets
44
45 // Muons
46 PromptFinalState bare_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
47 LeptonFinder all_dressed_mu(bare_mu, photons, 0.1, Cuts::abseta < 5.0);
48
49 // Electrons
50 PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
51 LeptonFinder all_dressed_el(bare_el, photons, 0.1, Cuts::abseta < 5.0);
52
53 //Jet forming
54 VetoedFinalState vfs(FinalState(Cuts::abseta < 5));
55 vfs.addVetoOnThisFinalState(all_dressed_el);
56 vfs.addVetoOnThisFinalState(all_dressed_mu);
57
58 FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
59 declare(jets, "Jets");
60
61 // Book auxiliary histograms
62 book(_h["MTWZ"], "_mTWZ", refData( 6, 1, 1));
63 book(_h["sumpt"], "_sumpT", refData( 8, 1, 1));
64 book(_h["dphiWZ"], "_dphiWZ", refData(10, 1, 1));
65 book(_h["mjj"], "_mjj", refData(14, 1, 1));
66 book(_h["dyjj"], "_dRapjj", refData(16, 1, 1));
67 book(_h["dphijj"], "_dPhijj", refData(18, 1, 1));
68 book(_d["Njets_VBS"], "_njets", refData<YODA::BinnedEstimate<string>>(12, 1, 1));
69 book(_d["Njets_gap"], "_gapJets", refData<YODA::BinnedEstimate<string>>(20, 1, 1));
70
71 // book output bar charts
72 book(_s["MTWZ"], 6, 1, 1);
73 book(_s["sumpt"], 8, 1, 1);
74 book(_s["dphiWZ"], 10, 1, 1);
75 book(_s["mjj"], 14, 1, 1);
76 book(_s["dyjj"], 16, 1, 1);
77 book(_s["dphijj"], 18, 1, 1);
78
79 }
80
81
82 void analyze(const Event& event) {
83
84 const Particles& dressedleptons = apply<LeptonFinder>(event, "LeptonFinder").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
88 int i, j, k;
89 double MassZ01 = 0., MassZ02 = 0., MassZ12 = 0.;
90 double MassW0 = 0., MassW1 = 0., MassW2 = 0.;
91 double WeightZ1, WeightZ2, WeightZ3;
92 double WeightW1, WeightW2, WeightW3;
93 double M1, M2, M3;
94 double WeightTotal1, WeightTotal2, WeightTotal3;
95
96 //---Fiducial PS: assign leptons to W and Z bosons using Resonant shape algorithm
97 if (dressedleptons.size() < 3 || neutrinos.size() < 1) vetoEvent;
98
99 int icomb=0;
100 // try Z pair of leptons 01
101 if ( (dressedleptons[0].pid() ==-(dressedleptons[1].pid())) && (dressedleptons[2].pid()*neutrinos[0].pid()< 0) && (dressedleptons[2].abspid()==neutrinos[0].abspid()-1)) {
102 MassZ01 = (dressedleptons[0].momentum() + dressedleptons[1].momentum()).mass();
103 MassW2 = (dressedleptons[2].momentum() + neutrinos[0].momentum()).mass();
104 icomb = 1;
105 }
106
107 // try Z pair of leptons 02
108 if ( (dressedleptons[0].pid()==-(dressedleptons[2].pid())) && (dressedleptons[1].pid()*neutrinos[0].pid()< 0) && (dressedleptons[1].abspid()==neutrinos[0].abspid()-1)) {
109 MassZ02 = (dressedleptons[0].momentum() + dressedleptons[2].momentum()).mass();
110 MassW1 = (dressedleptons[1].momentum() + neutrinos[0].momentum()).mass();
111 icomb = 2;
112 }
113 // try Z pair of leptons 12
114 if ( (dressedleptons[1].pid()==-(dressedleptons[2].pid())) && (dressedleptons[0].pid()*neutrinos[0].pid()< 0) && (dressedleptons[0].abspid()==neutrinos[0].abspid()-1)) {
115 MassZ12 = (dressedleptons[1].momentum() + dressedleptons[2].momentum()).mass();
116 MassW0 = (dressedleptons[0].momentum() + neutrinos[0].momentum()).mass();
117 icomb = 3;
118 }
119
120 if (icomb<=0) vetoEvent;
121
122
123 WeightZ1 = 1/(pow(MassZ01*MassZ01 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
124 WeightW1 = 1/(pow(MassW2*MassW2 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
125 WeightTotal1 = WeightZ1*WeightW1;
126 M1 = -1*WeightTotal1;
127
128 WeightZ2 = 1/(pow(MassZ02*MassZ02- MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
129 WeightW2 = 1/(pow(MassW1*MassW1- MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
130 WeightTotal2 = WeightZ2*WeightW2;
131 M2 = -1*WeightTotal2;
132
133 WeightZ3 = 1/(pow(MassZ12*MassZ12 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
134 WeightW3 = 1/(pow(MassW0*MassW0 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
135 WeightTotal3 = WeightZ3*WeightW3;
136 M3 = -1*WeightTotal3;
137
138 if( (M1 < M2 && M1 < M3) || (MassZ01 != 0 && MassW2 != 0 && MassZ02 == 0 && MassZ12 == 0) ) {
139 i = 0; j = 1; k = 2;
140 }
141 if((M2 < M1 && M2 < M3) || (MassZ02 != 0 && MassW1 != 0 && MassZ01 == 0 && MassZ12 == 0) ) {
142 i = 0; j = 2; k = 1;
143 }
144 if((M3 < M1 && M3 < M2) || (MassZ12 != 0 && MassW0 != 0 && MassZ01 == 0 && MassZ02 == 0) ) {
145 i = 1; j = 2; k = 0;
146 }
147
148 FourMomentum Zlepton1 = dressedleptons[i].momentum();
149 FourMomentum Zlepton2 = dressedleptons[j].momentum();
150 FourMomentum Wlepton = dressedleptons[k].momentum();
151 FourMomentum Zboson = dressedleptons[i].momentum()+dressedleptons[j].momentum();
152 FourMomentum Wboson = dressedleptons[k].momentum()+neutrinos[0].momentum();
153
154 double cosLepNeut;
155 double Wboson_mT = 0.;
156 double norm = Wlepton.pT() * neutrinos[0].pt();
157 if(norm != 0 ) {
158 cosLepNeut = ( Wlepton.px()*neutrinos[0].px() + Wlepton.py()*neutrinos[0].py() )/norm ;
159 if (1-cosLepNeut >= 0. ) Wboson_mT = sqrt( 2 * Wlepton.pT() * neutrinos[0].pt() * (1-cosLepNeut ) );
160 }
161
162 //---- CUTS (based on Table 1 WZ: 36.1 fb-1)----//
163 if (Wlepton.pT() <= 20*GeV || Zlepton1.pT() <= 15*GeV || Zlepton2.pT() <= 15*GeV) vetoEvent;
164 if (Wlepton.abseta() >= 2.5 || Zlepton1.abseta() >= 2.5 || Zlepton2.abseta() >= 2.5) vetoEvent;
165 if (fabs(Zboson.mass()/GeV - MZ_PDG) >= 10.) vetoEvent;
166 if (Wboson_mT <= 30*GeV) vetoEvent;
167 if (deltaR(Zlepton1, Zlepton2) <= 0.2) vetoEvent;
168 if (deltaR(Zlepton1, Wlepton) <= 0.3) vetoEvent;
169 if (deltaR(Zlepton2, Wlepton) <= 0.3) vetoEvent;
170
171 double WZ_pt = (Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt() + neutrinos[0].pt())/GeV;
172 double WZ_px = (Zlepton1.px() + Zlepton2.px() + Wlepton.px() + neutrinos[0].px())/GeV;
173 double WZ_py = (Zlepton1.py() + Zlepton2.py() + Wlepton.py() + neutrinos[0].py())/GeV;
174 double mTWZ = sqrt( pow(WZ_pt, 2) - ( pow(WZ_px, 2) + pow(WZ_py,2) ) );
175 double sumptleptons = (Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt())/GeV;
176 double dPhiWZTruth = acos(cos(Zboson.phi()-Wboson.phi()));
177
178
179
180 //---- Jet CUTS----//
181 idiscard(jets, [&](const Jet& j) {
182 return deltaR(j, Zlepton1) < 0.3 || deltaR(j, Zlepton2) < 0.3 || deltaR(j, Wlepton) < 0.3;
183 });
184 if (jets.size() < 2) vetoEvent;
185 if (jets[0].pT() < 40*GeV) vetoEvent;
186
187 // Selection of the second jet as the second highest pT jet and in opposite hemisphere with the fisrt jet
188 FourMomentum jet_lead = jets[0].mom();
189 FourMomentum jet_sublead;
190 bool foundVBSJetPair = false;
191 for (const Jet& jet : jets) {
192 if(jet.pT() > 40*GeV && jet.eta()*jets[0].eta() < 0.) {
193 jet_sublead = jet.mom();
194 foundVBSJetPair = true;
195 break;
196 }
197 }
198 if (!foundVBSJetPair) vetoEvent;
199
200 const double mJJ = (jet_lead + jet_sublead).mass()/GeV;
201 const double dphi_jj = acos(cos(jet_lead.phi() - jet_sublead.phi()));
202 const double dyjj = fabs(jet_lead.rap() - jet_sublead.rap());
203
204 //Plots in the SR
205 if (mJJ < 500*GeV) vetoEvent;
206
207 const size_t njets40 = select(jets, Cuts::pT > 40*GeV).size();
208 fillDiscrete("Njets_VBS", njets40, 5);
209
210 const double y_min = std::min(jet_lead.rap(), jet_sublead.rap());
211 const double y_max = std::max(jet_lead.rap(), jet_sublead.rap());
212 const size_t njetsGap = count(jets, [&](const Jet& j) {
213 return (j.rap() > y_min && j.rap() < y_max);
214 });
215 fillDiscrete("Njets_gap", njetsGap, 3);
216
217 fillWithOverflow("MTWZ", mTWZ, 551);
218 fillWithOverflow("sumpt", sumptleptons, 501);
219 fillWithOverflow("mjj", mJJ, 2001);
220
221 _h["dphiWZ"]->fill(dPhiWZTruth);
222 _h["dyjj"]->fill(dyjj);
223 _h["dphijj"]->fill(dphi_jj);
224
225 }
226
227 void fillWithOverflow(const string& tag, const double value, const double overflow) {
228 _h[tag]->fill(value < overflow? value : overflow);
229 }
230
231 void fillDiscrete(const string& tag, const size_t value, const size_t overflow) {
232 string edge = "$\\geq" + std::to_string(overflow) + "$";
233 if (value < overflow) edge = std::to_string(value);
234 _d[tag]->fill(edge);
235 }
236
237 void finalize() {
238
239 scale(_h, crossSectionPerEvent() / femtobarn);
240 scale(_d, crossSectionPerEvent() / femtobarn);
241 // unfortunately, no differential cross-sections were measured in this analysis
242 for (auto& item : _h) barchart(item.second, _s[item.first]);
243
244 }
245
246
247 /// @}
248
249 private:
250
251
252 /// @name Histograms
253 /// @{
254 map<string, Histo1DPtr> _h;
255 map<string, BinnedHistoPtr<string>> _d;
256 map<string, Estimate1DPtr> _s;
257 /// @}
258
259 const double MZ_PDG = 91.1876;
260 const double MW_PDG = 80.385;
261 const double GammaZ_PDG = 2.4952;
262 const double GammaW_PDG = 2.085;
263
264 };
265
266
267 RIVET_DECLARE_PLUGIN(ATLAS_2018_I1711223);
268
269}
|